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A surrogate model methodology is described for predicting in real time the resid- 
ual strength of flight structures with discrete-source damage. Starting with design 
of experiment, an artificial neural network is developed that takes as input discrete- 
source damage parameters and outputs a prediction of the structural residual strength. 
Target residual strength values used to train the artificial neural network are derived 
from 3D finite element-based fracture simulations. Two ductile fracture simulations are 
presented to show that crack growth and residual strength are determined more accu- 
rately in discrete-source damage cases by using an elastic-plastic fracture framework 
rather than a linear-elastic fracture mechanics-based method. Improving accuracy of 
the residual strength training data would, in turn, improve accuracy of the surrogate 
model. When combined, the surrogate model methodology and high fidelity fracture 
simulation framework provide useful tools for adaptive flight technology. 
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E = elastic modulus (GPa) 
u = Poisson’s ratio (mm/mm) 

(jy = yield stress (MPa) 

STRI65 = quadratic triangular shell element in ABAQUS [1] 

SSR = quadratic reduced-integration shell element in ABAQUS [1] 

C3D10 = quadratic tetrahedral elements in ABAQUS [1] 

C3D15 = quadratic wedge element in ABAQUS [1] 

C3D20{R) = quadratic brick element (reduced-integration) in ABAQUS [1] 
a = crack length (cm) 

n = number of cracks in discrete-source damage 

0 = orientation of discrete-source damage, angle between positive x axis and nearest crack 
dx = distance from middle stiffener to center of discrete-source damage (cm) 

Ej , Kii, Kill = mode I, II, and III plane strain stress intensity factors (MPa^/m) 

Kic^Kjic = plane strain fracture toughness for modes I and II (MPa^/m) 

Pmax = damage-dependent allowable traction, residual strength (MPa) 

P = applied traction (MPa) 

MSE = mean squared error as defined by Eq. (2) 

Cy = correlation coefficient as defined by Eq. (3) 

CTD = magnitude of relative displacement between upper and lower fracture surfaces (mm) 
CTDcrit = critical value of CTD (mm) 

CTD I, CTDii, CTD III = opening, in-plane sliding, out-of-plane shearing components of CTD 
(mm) 

d = fixed characteristic distance behind crack front where CTD is monitored (mm) 
da = crack extension (mm) 

(n) = script to denote mesh at crack increment 

(n + 1) = script to denote mesh at (n -h 1)^^ crack increment 
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I. Introduction 


Resilient aircraft control involves adaptive responses to off-nominal flight conditions, including 
the incurrence of structural discrete-source damage during flight. Discrete- source damage is typically 
manifested as a result of a structural impact event, including hail- and birdstrike. In 2003, an Airbus 
A300 operated by DHL was struck by a surface-to-air missile after takeoff from Baghdad, Iraq, 
causing discrete-source damage to crucial control surfaces of the left wing [2]. In 2008, a Boeing 
747-438 operated by Qantas Airways incurred in-flight structural damage to the fuselage and right 
wing leading edge following the failure of an onboard oxygen cylinder [3]. Although the aircraft 
landed safely in both cases, these examples motivate a need for more resilient, adaptive control 
system responses. 

In these types of cases, problems associated with in-flight discrete-source damage, for example 
inability to sustain original design loads, can be exacerbated by crack propagation from damaged 
regions. To avoid unstable crack propagation, load levels must be maintained below a reduced 
load-carrying capacity, or residual strength^ of damaged flight structures. Adaptive control system 
responses might include automatic adjustment of certain flight parameters (e.g. velocity, maximum 
acceleration) to accommodate structural residual strength. This accommodation implies that accu- 
rate residual strength predictions of flight structures with complex damage configurations be made 
in real time, during flight] this capability currently does not exist for commercial aviation. 

Challenges to developing an adaptive response technology include accurately predicting residual 
strength of discrete-source damaged structures both offline (i.e. during control system design) and 
online (i.e. in real time onboard the aircraft). In the offline context, researchers have developed 
various tools for determining residual strength of thin, damaged metallic structures using elastic- 
plastic fracture mechanics (EPFM)-based numerical methods. For example, two common finite 
element (FE) modeling techniques involve nodal release and adaptive remeshing. Both techniques 
represent cracks geometrically [4]. The former, however, prescribes possible crack trajectories, which 
introduces inherent mesh dependencies into fracture simulations and limits generality of crack path 
predictions. Nodal release techniques have been used in 2D [5-12] and in 3D [13-15] for studying 
crack growth parameters and predicting residual strength of structures where the crack path was 
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known a priori and where mesh refinement along the crack path sufficiently characterized growth 
increments. Adaptive remeshing techniques avoid such mesh dependencies and enable simulation 
of arbitrary crack propagation using evolutionary models or criteria [16-20]. Adaptive remeshing 
techniques have been implemented in both 2D [21] and in 3D [22, 23]. Of the described techniques, 
3D, adaptively remeshed, elastic-plastic tearing simulations provide the most general prediction 
capabilities for crack growth and residual strength. 

It is infeasible to perform a rigorous and computationally intensive crack growth simulation 
within the possible short time span following a discrete-source damage event. Thus, an approxima- 
tion, or surrogate models is needed for making online predictions of residual strength. Queipo et al 
provided a complete description of surrogate modeling development and optimization [24]. With 
regard to surrogate construction, they described both parametric (e.g. polynomial regression and 
Kriging) and nonparametric (e.g. radial basis functions) approaches. In nonparametric approaches, 
a global functional form relating system input to system response is not assumed. 

Artificial neural networks (NNs) are a nonparametric surrogate modeling approach and are 
trained to infer a nonlinear mapping from system input to system response, or output. The reader is 
referred to [25] for an extensive methodology overview of the most widely used types of NN. Different 
types of NNs have been applied extensively for damage detection [26-32] and, to a much lesser extent, 
for damage assessment. Queues et al. employed a NN methodology to predict fracture indicators 
(e.g. density of fractures) in naturally fractured rock reservoirs as a function of various geological 
and geophysical data [33]. Pidaparti et al. employed a NN to predict residual strength and corrosion 
rate of aging aircraft panels with collinear multi-site damage by training with experimental results 
and validating with both experimental results and analytical solutions [34]. Recently, Mohanty et 
al. used a Gaussian process (GP) approach to predict fatigue crack growth in aluminum 2024-T351 
specimens by training two distinct models, one presented with experimental load parameters as 
input and another presented with piezoelectric sensor signals as input [35]. In that work, Mohanty 
et al. used observed fatigue crack lengths and growth rates as known output for training each model. 

Alternatively, NNs can be trained using results from numerical experiments, or simulations [36]. 
For example, Sankararaman et al. recently used linear-elastic fracture parameters computed from 
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FE analyses to train a GP model as part of a method to statistically infer equivalent initial flaw 
size in fatigue applications [37]. High-fldelity numerical simulations can provide training data when 
analytically- and experimentally-derived data is limited due either to a lack of generally applicable 
analytical solutions or to prohibitive costs of obtaining sufficient experimental data. 

The purpose of the work presented here is two-fold: (1) to illustrate a methodology for creating 
a surrogate model as a real-time residual strength prediction tool and (2) to describe and validate 
numerical tools for making accurate residual strength predictions offline using fully 3D, elastic- 
plastic, FE-based crack growth simulations. The high-fidelity, more computationally expensive tools 
described in (2) can provide training data that, when coupled with the surrogate model methodology 
described in (1), can be used in the design of adaptive response technology. 

Consistent with our two-fold purpose, this paper is divided into two primary sections. Section II 
illustrates the methodology for developing a surrogate model (in particular, a NN) that predicts 
residual strength as a function of discrete-source damage parameters. The methodology is illustrated 
using a relatively simple proof-of-concept example. The procedure for gathering training data 
is described in II A and II B. Because an implementation-ready NN is beyond the scope of this 
paper, training data for the proof-of-concept example relies on reduced-order residual strength 
approximations. After collecting training data, a simple NN is constructed in II C by optimizing 
certain performance parameters. Finally, a sensitivity study is conducted in II D to understand the 
effect of each damage parameter on predicted residual strength specifically for the proof-of-concept 
structure. 

Section III improves upon offline residual strength prediction tools used in Section II by simulat- 
ing 3D, elastic-plastic tearing. The tools provide more general crack growth simulation capabilities 
and can be used to generate accurate residual strength training data. Simulation results are vali- 
dated in III C for a mixed-mode I/II fracture test and for a relatively large, integrally-stiffened panel 
(ISP) that exhibits crack branching. 

Results and discussions from the NN proof-of-concept example and from the elastic-plastic tear- 
ing simulations are provided in each respective section. Section IV offers a summary and conclusions 
for the entirety of this work. 
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II. Neural Network Development and Methodology 


This section describes the development of a surrogate model for predicting residual strength of 
discrete-source damaged aircraft structures in real time. A global functional form is not assumed 
for the nonlinear relationship between residual strength and the damage parameters influencing it; 
thus, a nonparametric surrogate model is developed. In particular, a supervised NN is considered 
due to rapid prediction capabilities amenable to real-time applications. In Fig. 1, the upper dashed 
region shows the generalized procedure for developing a NN (surrogate model) that predicts residual 
strength as a function of parameterized discrete-source damage. The lower dashed region shows the 
functionality of the NN (surrogate model) in a real-time context. 



Fig. 1 Upper dashed box illustrates a general approach for developing a surrogate model to 
predict residual strength of damaged structures. Lower dashed box illustrates how the sur- 
rogate model would function onboard an aircraft for predicting residual strength of damaged 
structures in real time. 


The first step in this type of surrogate model development is typically referred to as design of 
experiment (DOE) [24] and involves obtaining data points that will be used to train and test the 
NN. The DOE should be based on the intended application of the NN. For example, if the NN is 
intended to provide residual strength predictions in terms of maximum allowable bending moment 
in a damaged aircraft wing, then the data points should be gathered using an appropriate wing 
structure with applied boundary conditions of interest. Each data point includes sampled input 
variable(s) and corresponding known system response(s), called target output. Once the NN has 
been trained to map given input to target output, it becomes a useful tool for predicting system 
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response when presented with new input that is within the training range but does not necessarily 
correspond to data points used for training. 

To illustrate the methodology, a simple NN is developed using a representative wing structure 
and reduced-order (linear-elastic) approximations for predicting residual strength. The representa- 
tive wing structure is a 61.0 x 91.4 cm^ integrally-stiffened panel (ISP) with three blade stiffeners 
each 5.1 cm in height, as shown in Fig. 2. The ISP skin and stiffeners are 2.3 mm thick. The panel 
is modeled as linear-elastic with = 71.0 GPa and v = 0.33, similar to values for a 2XXX series 
lower wing skin aluminum alloy (AA). 


P 



Gx I I I l5.1 


(thickness of skin and stiffeners = 0.23 ) 


Fig. 2 Schematic of ISP model with dimensions similar to those in [38]. Plan view (top) 
and cross-section showing integral blade stiffeners (bottom). A damage-containing region is 
modeled using 3D solid elements (enclosed in shell-solid boundary) while remainder of panel 
is modeled with shell elements. All dimensions in cm. 


Multiple FE models of the uncracked panel are constructed using ABAQUS(R) [1]. A shell-solid 
modeling technique is employed, where each panel is modeled using 3D solid elements in a region 
that will contain damage and shell elements elsewhere, as depicted schematically in Fig. 2. In this 
way, 3D constraint is inherently captured along crack fronts using fully 3D solid elements, while 
shell elements help maintain a level of computational efficiency without losing capability to capture 
out-of-plane deformation and possible buckling. The shell and solid element regions are joined 
using a coupling constraint, whereby resultant forces and moments acting at shell edge nodes on the 
shell-solid boundary are distributed as forces acting at nodes located in a region of influence on the 
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solid surface of the shell-solid boundary. A mesh refinement study is carried out to ensure adequate 
discretization of the panel models. Uncracked panels are modeled using approximately 50 S TRIO 5, 
2000 S8R, and between 1800 and 17,300 C3D20R elements, depending on the size of the damaged 
region. Boundary conditions for the ISP models are defined to emulate tensile loading conditions 
for a region of the lower wing surface and are shown schematically in Fig. 2. 

A supplementary study was carried out to determine shell-solid boundary effects on nearby 
crack fronts in order to minimize the size of the solid region without affecting stress intensity factors 
(SIFs) computed along nearby crack fronts. Maintaining fracture parameter accuracy is especially 
important since fracture parameters are used to predict structural residual strength (described in 
II B), which is in turn used to train the NN (described in II C). The supplementary study considered 
a 61.0 X 91.4 cm^ unstiffened panel of the same (linear-elastic) material and thickness as the ISP 
described above. The panel had a single, 12.7 cm long, centrally-located through-crack oriented in 
the X direction (normal to applied tensile load). Both tensile and bending conditions were considered 
in the study. The panel was modeled entirely with shell elements except for a region containing the 
crack, which was modeled with 3D solid elements. All model parameters remained constant while 
varying the size (both in-plane dimensions) of the square-shaped solid region, and therefore the 
distance from the shell-solid boundary to the crack front. The size of the solid region was initially 
slightly larger than the crack length and was increased until computed SIFs converged. Results from 
the supplementary study indicated that for a static, linear-elastic crack, the distance from shell-solid 
boundary to nearest crack front should be no less than 25% of the crack length. This ensures that 
the shell-solid boundary has negligible effect on computed SIFs. The same rule-of-thumb is applied 
to the example ISP models described in the NN study. 

The following sections describe the generally applicable methodology for developing a NN as a 
real-time residual strength prediction tool. 

A. Input Variables: Discrete-source Damage Parameters 

Discrete-source damage in this work is represented by a symmetric, star-shaped, array of equi- 
length cracks, as depicted in Fig. 3(b). This representation of discrete-source damage is motivated 
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by observations of petaling caused by penetration damage to thin metallic structures, see Fig. 3(a). 


If all of the cracks in the star-shaped array of Fig. 3(b) separate under load (i.e. there are no crack 


closure effects), then the cracked region transfers no load and effectively represents a circular hole 


with petaling edges, similar to that shown in Fig. 3(a). The damage representation is parameterized 


by the four variables n, a, and d, which are postulated to influence residual strength of the ISP. 



Reprinted with permission of 
ASM International®. 

All rights reserved. 
www.asminternational.orq 


a) 



Fig. 3 (a) Petaling on the reverse side of a metallic sheet subject to explosive, discrete-source 
damage [39]. (b) Schematic showing the representation and parameterization of discrete- 

source in the NN example described in this work. 


The sample space of damage configurations is defined by a range of values for each parameter. 
Ranges can be specified based on accident reports, photographic evidence, potential structural 
threats, design specifications, and so forth. Inherently, the NN predictions are valid only for input 
parameter values within the range of training data. Thus, it is necessary to define the sample space 
based on the particular NN application. In the example NN, ranges for each damage parameter 
are limited to some extent by the ISP geometry. Each range is given in Table 1. The parameter 
n, takes integer values ranging from two to six. The range of 0 depends on n due to the definition 
of orientation and the symmetry of the star-shaped configuration. The range of a is defined in 
terms of ISP bay width, from 1/8 * hay width to 1/4 * bay width. Due to symmetry of the ISP 
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model, the parameter dx ranges from 10.2 cm (damage centered in mid-bay) to 0 (damage centered 
at middle stiffener). If the damage is located such that the damage-containing, solid FE region 
overlaps anywhere with the middle stiffener, the stiffener is assumed to be severed in the damaged 
region and is modeled explicitly as such. 


Table 1 Range of values associated with each damage parameter in the example NN. 


Damage Parameter Range 


n 

2-6 

n = 2 (deg) 

0-90 

n = 3 (deg) 

0-60 

n = 4 (deg) 

0-45 

n = 5 (deg) 

0-36 

n = 6 (deg) 

0-30 

a (cm) 

1.27-5.08 

dx (cm) 

0-10.2 


The damage parameter space is sampled to obtain damage configurations, each expressed as 
a combination of input parameters (n^O^a^dx). The space of variables can be sampled using a 
number of different sampling methods, including random, stratified, and Latin Hypercube [40]. 
Latin Hypercube Sampling (LHS) is a type of stratified sampling method that guarantees each 
partition, or stratum, of input variable space is sampled, though not necessarily uniformly. In this 
work, LHS is performed five times for each of the variables (6>, a, and dx). Each of the five LHS runs 
corresponds to a different value of n (two, ..., six cracks) and requires the number of partitions to 
be specified. The MATLAB@ implementation for LHS is used here [41], where output is provided 
in the range from zero to one. Each sample value is then scaled to the respective parameter range 
according to Table 1. 

Table 2 shows all damage configurations (26 in total) that are modeled in the ISP NN example, 
where each configuration is expressed in terms of sampled input parameters. For each damage 
configuration, the x and y dimensions of the square, damage-containing, solid FE region are provided 
in the sixth column. The x and y dimensions are each 25% larger than the diameter of the star- 
shaped damage (i.e. 1.25 * 2a), as suggested by the supplementary shell-solid boundary effect study 
described above. The last column specifies whether or not the solid, damaged region severs the 
middle stiffener. If so, the portion of the stiffener that intersects the solid model region is removed; 
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otherwise the stiffener remains intact. 


Table 2 Damage configurations modeled in the ISP NN example. Each damage configuration 
is assigned an alphanumeric identification with number corresponding to n. The sixth column 
provides x and y dimensions of the square region in the shell-solid ISP. 


Damage 

configuration 

a dx 0 

ID (cm) (cm) (deg) 

n Solid region Severs 

x^y dimensions (cm) stiffener? 

2A 

3.8 

7.1 

21.8 

2 

9.39 

NO 

2B 

4.2 

2.1 

87.8 

2 

10.48 

YES 

2C 

3.0 

3.6 

2.2 

2 

7.53 

YES 

2D 

1.4 

0.2 

35.7 

2 

3.49 

YES 

2E 

2.3 

9.9 

36.5 

2 

5.66 

NO 

2F 

4.5 

6.1 

44.9 

2 

11.25 

NO 

3A 

3.7 

8.3 

5.0 

3 

9.36 

NO 

3B 

1.5 

0.6 

25.2 

3 

3.72 

YES 

3C 

2.9 

9.6 

23.2 

3 

7.15 

NO 

3D 

4.0 

3.7 

32.9 

3 

9.88 

YES 

3E 

4.6 

1.7 

10.4 

3 

11.56 

YES 

3F 

2.0 

6.3 

38.6 

3 

4.92 

NO 

4A 

1.7 

6.5 

6.6 

4. 

4.15 

NO 

4B 

3.4 

1.7 

18.7 

4 

8.60 

YES 

4C 

4.9 

9.7 

25.1 

4 

12.3 

NO 

4D 

2.1 

4.4 

27.0 

4 

5.37 

NO 

4E 

3.0 

7.0 

14.9 

4 

7.52 

NO 

5A 

3.3 

8.2 

4.8 

5 

8.30 

NO 

5B 

2.2 

0.8 

6.9 

5 

5.43 

YES 

5C 

1.5 

5.4 

19.8 

5 

3.84 

NO 

5D 

3.2 

9.4 

22.6 

5 

7.91 

NO 

6A 

1.8 

8.3 

5.4 

6 

4.50 

NO 

6B 

3.7 

9.7 

26.5 

6 

9.21 

NO 

6C 

4.9 

6.0 

12.2 

6 

12.20 

NO 

6D 

4.2 

2.0 

18.4 

6 

10.61 

YES 

6E 

3.0 

3.5 

21.9 

6 

7.54 

YES 


B. Target Output: Residual Strength from Numerical Fracture Simulations 

For each input damage configuration, a numerical fracture simulation is employed to determine 
residual strength, which provides target output used to train and test the NN. FRANC3D\NG [42] 
is used to insert each parameterized star-shaped crack configuration into the solid FE region of each 
panel. An ABAQUS@ contact algorithm is employed to prevent crack surfaces from overlapping 
during the applied loading. Contact properties are defined as frictionless in the tangential direction 
with “hard” pressure-overclosure behavior normal to the contacting crack surfaces, which minimizes 
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interpenetration. The FE models are then analyzed using ABAQUS@, and FE analysis results are 
post-processed to determine residual strength. 

For the sake of illustrating the NN methodology, two simplifying assumptions are made here 
to predict residual strength of the ISPs. First, the ISPs remain linear-elastic and can be analyzed 
using LEFM parameters (SIFs). Second, the residual strength can be predicted for a static crack 
configuration (i.e. crack growth is not modeled in this example). 

In the ISP NN example, the LEFM approximation of residual strength is based on mixed-mode 
I/II fracture criteria [16, 18] to account for local mode mixity (in-plane) of angled cracks in the 
star-shaped damage array. In [43], Broek describes a practical mixed-mode I/II failure envelope, 
approximated by the equation of an ellipse: 


{KilKi.f + {KnlKucf = 1 - ( 1 ) 

For the AA 2XXX series material in the ISP example, = 32 MPay^, and Kn^ is assumed to 
be 10% less than Kic after results from the strain energy density criterion presented by Sih [18]. 

Using this mixed-mode LEFM-based approximation, residual strength is defined here as the 
applied traction load. Fig. 2, that first causes unstable crack growth for any point along any crack 
front of the star-shaped damage configuration. In other words, as soon as one point along one crack 
front reaches a critical combination {Ki,Kn)c on the elliptical failure envelope, the entire panel is 
assumed to fail. The method for determining the residual strength for each damaged panel is shown 
in Fig. 4 and proceeds as follows: (1) analyze the ISP FE model with P; (2) compute Kj and Ku 
at each node along each crack front using FRANC3D\NG; (3) for each crack front node, find the 
intersection point {Kj, Ku)c of the elliptical failure envelope with a straight line from the origin to 
the computed {Ki,Ku) and subsequently find the linear scaling factor. A, that maps {Ki,Kn) to 
{Kj,Kjj)c] (4) of all the computed scaling factors, select the most critical, Ac; (5) calculate Pmax 
as P scaled by Ac. 

To ensure that nonlinearity due to crack face contact does not invalidate the linear load scaling 
approach described above, each of the damaged ISPs is reanalyzed with the respective scaled load. 
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crack face contact 



Fig. 4 LEFM-based procedure for approximating residual strength of the damaged ISPs in 
the NN example. 


i.e. the approximated residual strength. In all cases, {Kj^Ku) = {Kj^Kjj)c at the predicted crack 
front failure point, indicating that the scaled loads indeed correspond to failure loads according to 
the LEFM-based failure criterion assumed for this example problem. Values of Pmax provide the 
target outputs used to train the NN. 


C. Neural Network Construction 

The inputs (sampled damage parameters) and target outputs (residual strength predictions from 
numerical fracture simulations) are used to train and test a NN. For the ISP example, a feedforward 
NN with a backpropagation learning rule [25, 44], which is a commonly used type of supervised 
NN, is constructed using MATLAB@ [41]. The NN consists of a single hidden layer mapping the 
four-parameter input vectors (n, a, d^) to the single- valued outputs (Pmax)- The reader is referred 
to [44] for a general discussion on NNs and details regarding specific implementation of the transfer 
functions and training algorithm described next. A tan-sigmoid transfer function is employed to 
map the weighted inputs plus bias to the interval (-1,1)* ^ linear transfer function proportionally 
maps the weighted output plus bias from the hidden layer to the output layer. Data presented 
to the NN is divided into three sets — training, validation, and test. Weights and biases of the 
NN are adjusted at each iteration, or epoch, using the training set and a Levenberg-Marquardt 
optimization algorithm, as described in [45]. The algorithm seeks to improve performance of the 
NN by minimizing error between the NN outputs and the target outputs. Weights and biases from 
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training at any epoch are then used to check performance of the NN using the validation and test 
sets. The validation set prevents overtraining of the NN by ceasing training if performance degrades 
over a certain number of successive epochs. The test set is not used for training but is used to test 
NN accuracy following the current training epoch. The NN performance metric used here is the 
MSE, calculated as: 

= ( 2 ) 

where the superscript (i) corresponds to the training, validation, or test set, Q is the number of 
data points in the respective set, tk is target output for the input, and Pk is output predicted 
by the NN for the same input. 

The NN can be optimized by adjusting any number of parameters, including transfer functions 
between layers, number of hidden layers, various performance metrics, and so forth. In the ISP 
example, the NN is optimized by varying the number of neurons in the hidden layer (4,5,6) and by 
increasing size of the training set from 60%, to 70%, to 80% of the available data (with the balance 
equally divided between validation and test sets). Further, the performance metrics are optimized 
by minimizing MSE for the training and testing sets and by specifying that the correlation coefficient 
between NN output and targets should be at least 0.95 over the entire data set. 

D. Parametric Sensitivity Studies 

The trained NN is then employed to conduct parametric sensitivity studies, whereby sensitivity 
of residual strength to each postulated damage parameter is gauged. The sensitivity studies are 
carried out for configuration 4E, Table 2, as it represents an average damage configuration in terms 
of n, a, and dx as compared to the other configurations. 

A sensitivity study is carried out for each of the four damage parameters. In each study, 
three damage parameters of configuration 4E are held constant while one is varied. The variable 
parameter in each study takes values in the range of the corresponding variable on which the NN 
was trained. For example, the longest a considered in the sensitivity study is no longer than the 
longest a used to train the NN, which is a = 4.9 cm in the ISP example (damage configurations 
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4C and 6C). Further, the variable parameter takes values that are equally incremented within the 
respective range. Results from the study are presented in the following subsection. 

E. Results and Discussion from Neural Network Example 

Table 3 shows the approximated residual strengths of all damaged ISPs based on numerical 
fracture analyses and LEFM assumptions outlined in IIB. The table is sorted in order of increasing 
residual strength, and the corresponding damage parameters are provided to help draw preliminary 
conclusions. One immediate observation is that panels with severed stiffeners have lower (~ 50 — 
80 MPa) residual strengths, as expected. The single exception is damage configuration 2B. Though 
it severs the stiffener, configuration 2B is less critical than all other stiffener-severing cases and some 
intact-stiffener cases because it is an n = 2 configuration (straight crack) aligned with the loading 
direction. Overall, the correlation between severed stiffener and reduced residual strength highlights 
the effect of the load carrying stiffener on crack criticality. 

The ISPs with the lowest computed residual strength (configuration 3D) and highest computed 
residual strength (configuration 5C) are presented in Fig. 5. Each ISP is depicted with its respective 
residual strength, or failure load, applied. The predicted point of first-failure lies along the crack 
front indicated. Configuration 5C has more cracks and is 62.5% smaller than configuration 3D, 
though it is not the smallest of all configurations. More importantly, configuration 5C leaves the 
stiffener intact while configuration 3D results in a severed stiffener. For configuration 3D, the crack 
front that lies within the severed region and near the geometric discontinuity of the stiffener junction 
is subjected to higher stresses and is predicted to be critical. 

The optimal NN consists of four neurons in a single hidden layer with 80% of available data (i.e. 
twenty damage configurations) allocated to training. NN performance metric (MSE) as a function 
of training epochs is plotted in Fig. 6 for training, validation, and test sets. The NN is best trained 
at epoch 141, beyond which the MSE in the validation set continually increases and overtraining is 
said to occur. At this epoch, MSE of the three sets are = 0.001, = 0.87, and 

_ Q 3 Q Weights and biases connecting the input layer (damage parameters) to the hidden 
layer and the hidden layer to the output (residual strength) at training epoch 141 are presented in 
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Table 3 LEFM-based residual strength approximations for all damage configurations consid- 
ered in the ISP example, sorted by increasing residual strength. 


Damage 

configuration 

a dx Severs Pmax 

ID (cm) (cm) stiffener? (MPa) 

3D 

4.0 

3.7 

YES 

54.2 

6D 

4.2 

2.0 

YES 

56.8 

3E 

4.6 

1.7 

YES 

56.8 

4B 

3.4 

1.7 

YES 

58.1 

6E 

3.0 

3.5 

YES 

63.5 

2C 

3.0 

3.6 

YES 

67.6 

5B 

2.2 

0.8 

YES 

71.2 

3B 

1.5 

0.6 

YES 

73.8 

2D 

1.4 

0.2 

YES 

80.7 

6C 

4.9 

6.0 

NO 

82.7 

4C 

4.9 

9.7 

NO 

83.9 

2A 

3.8 

7.1 

NO 

89.1 

3A 

3.7 

8.3 

NO 

90.3 

2B 

4.2 

2.1 

YES 

93.8 

5A 

3.3 

8.2 

NO 

94.3 

4E 

3.0 

7.0 

NO 

98.4 

5D 

3.2 

9.4 

NO 

101.5 

2F 

4.5 

6.1 

NO 

102.1 

6B 

3.7 

9.7 

NO 

109.2 

3C 

2.9 

9.6 

NO 

110.2 

4A 

1.7 

6.5 

NO 

124.9 

3F 

2.0 

6.3 

NO 

128.2 

2E 

2.3 

9.9 

NO 

129.5 

6A 

1.8 

8.3 

NO 

130.7 

4D 

2.1 

4.4 

NO 

132.8 

5C 

1.5 

5.4 

NO 

146.7 


Tables 4 and 5. 

Considering the entire set of damage configurations, the NN predictions correlate well with 
the target outputs at epoch 141, as depicted in Fig. 6. Despite the good overall correlation and 
small MSE for the training set, the MSE in the validation and testing sets (which include only 
three damage configurations each) may be too large for actual implementation onboard an aircraft, 
depending on design specifications. It is suspected that adding more samples to the entire set of 
damaged ISPs would further reduce these errors in the NN. 

The infiuence of each damage parameter on predicted residual strength can be visualized graph- 
ically by plotting predicted residual strength as a function of each damage parameter (see Fig. 7). 
Sensitivity can be quantified by a number of different metrics, many of which yield comparable 
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Fig. 5 Two different damaged panels (ID 5C and ID 3D) shown with respective Pmax applied. 
Panels represent damage configurations that are least critical (a) and most critical (b) of 
all configurations considered. Predicted failure point lies along the indicated crack front. 
Deformation is scaled by factor of 10. FE mesh is not shown for better contour visualization. 


Table 4 NN weights and biases used to map input layer to hidden layer for the optimized NN 
at training epoch 141. 



Input 

parameter 

Hidden layer neuron 
12 3 4 


n 

-0.43 -2.98 -1.11 1.83 

Weights 

a 

0.69 -9.17 1.61 -2.27 


X 

3.58 -0.15 1.50 -8.74 


0 

-1.53 2.87 -2.96 5.58 

Biases 


-0.42 0.45 1.37 2.94 


results [46]. Here, sensitivity is quantified by the c^, expressed as a percentage: 


Cy 


where a = 


lElliPrr 


-Prr^ 


N -1 


(3) 


For any given sensitivity subset, i corresponds to the sample configuration, Pmax is the average 
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Fig. 6 (a) NN performance as a function of training epochs for the optimized NN; overtraining 
occurs after epoch 141. (b) Correlation between predicted and target residual strength values 
considering all damage configurations. 


Table 5 NN weights and bias used to map hidden layer to output for the optimized NN at 
training epoch 141. 



Hidden layer neuron 

Output 


1 

0.78 

Weights 

2 

0.25 


3 

-1.70 


4 

1.05 

Bias 


1.29 


residual strength of the subset, and N is the total number of damage configurations in the respective 
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subset. Sensitivities to each damage parameter are calculated as = 24.8%, = 16.6%, 

= 6.0%, and = 1.2%. 

Orientation and number of cracks are found to have relatively minor influences on predicted 
residual strength, which is apparent both by their sensitivity metrics and by the plots (b) and (d) 
of Fig. 7. Crack length, on the other hand, has a more significant influence and causes a reduction 
in predicted residual strength as crack length increases, which is expected. What is unexpected, 
however, is the step-like behavior depicted in Fig. 7(c). This behavior is caused by binary modeling 
of the stiffener (explicitly modeling as severed or intact), a feature that is inherently implicit in 
both crack size and location. The stiffener effect is also apparent in Fig. 7(a) of damage location 
sensitivity. Predicted residual strength is lowest (and relatively insensitive to damage location) if 
the damage is located such that it severs the stiffener. As the damage location moves away from the 
stiffener and is no longer severing it, there is a linear increase in residual strength until the damage 
is located within the middle quarter of the bay. 

In general, the importance of this kind of sensitivity study is (1) to gain a better intuition of how 
and why certain damage characteristics influence residual strength and (2) to potentially decrease 
the dimensionality of the NN by neglecting parameters deemed insignificant. 

III. 3D Elastic-plastic Fracture Simulations for Improved Neural Network Training 

For discrete-source damage cases involving significant ductile tearing, a generally applicable 3D 
EPFM framework may be used for improving residual strength training data. An elastic-plastic 
crack growth simulation procedure, as implemented in this section, is illustrated in Fig. 8 and 
proceeds as follows: 

1. define an uncracked FE model and boundary conditions; 

2. extract a sub-region of the mesh for crack insertion, remeshing, and reconnection with the 
global mesh; 

3. map previous deformation and material state onto the remeshed model; 

4. perform nonlinear FE analysis and monitor the crack growth criterion; 
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Fig. 7 Sensitivity of predicted residual strength to each damage parameter. 


5. once criterion is satisfied, stop the current FE analysis to update crack configuration, remesh 
sub-region, and reconnect sub-region mesh with global mesh; 

6. repeat from step 3 until critical crack length is achieved or until residual strength is attained. 

The simulation procedure allows for prediction of curvilinear crack paths and arbitrary crack 
front evolution. The EPFM framework was overviewed in [47] and is described here for completeness. 
Scripts used for implementation are found in the Appendix. 


A. Nonlinear Fracture Parameter: Crack-tip Displacement 

In elastic-plastic tearing simulations, especially of thin metallic structures, crack growth should 
be characterized by an appropriate nonlinear parameter. One such parameter arises from correlation 
between crack growth and a critical amount of opening or displacement behind the crack tip (see 
[48] for details). A criterion based on this parameter, which is called the crack-tip displacement 
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*lnsert crack(s) or predict 
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*Map material state 
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Fig. 8 Elastic-plastic crack growth simulation algorithm using FRANC3D\NG. Contributions 
from this work include evaluation of crack-tip displacement (CTD) criterion during nonlin- 
ear FE analysis and implementation of material state mapping algorithm following adaptive 
remeshing. 


[CTD) or sometimes referred to as the generalized crack (tip) opening displacement [22, 49], is 
implemented here. Notably in simulation, once a critical value, CTDcrit , has been determined for 
a specific material and thickness through a calibration procedure, the same CTDcrit is applicable 
over a range of structural configurations comprising the same material and thickness under similar 
loading. In the EPFM simulations here, CTD is computed as: 

CTD = ^CTDi^ + CTDii^ + CTDiu^ (4a) 


CTD I = vi — V2 


(4b) 


CTD II = u\ — U2 


(4c) 


CTD III = wi~ W2, 


(4d) 


where v, and w correspond to displacements in the x, y, and 2 ) directions, respectively, and 


subscripts 1 and 2 denote the two points used to compute CTD. CTD is computed between two 


points that are initially coincident (one on each crack face) on the undeformed crack surface at a 


distance, d, behind a crack front node (i.e. in the direction normal to the crack front at the particular 
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crack front node). This is illustrated schematically in 2D in Fig. 9. In 3D, CTD values are computed 
behind multiple crack front nodes. The pair of initially coincident points where CTD is calculated 
is called a CTD point. Element shape functions are used to interpolate displacements {u^v^w) such 
that the CTD points need not correspond to nodal locations. 



— previous crack surface 
new crack surface 

o initially coincident points where CTD is computed 
Fig. 9 Simplified schematic of CTD implementation illustrated on crack profile. 

Crack growth occurs when CTD attains a critical value, CTDcritj within a specified tolerance. 
There are several ways to evaluate the CTD criterion when modeling a 3D crack front, including 
evaluation at a single CTD point either midway along the crack front or on the specimen’s free sur- 
face. Alternatively, the CTD criterion may be evaluated by comparing CTDcrit to an average CTD 
value calculated for multiple CTD points. Because CTDcrit is known to depend on 3D constraint 
at any point along a crack front, using a single CTDcrit to predict the advance of an entire crack 
front might not be valid for all cases. A more rigorous and computationally expensive evaluation 
technique would be to compare CTD at each CTD point to a constraint-dependent CTDcrit- While 
some work has been done to resolve a relationship between 3D constraint and CTDcrit [20, 50], a 
constraint-dependent fracture criterion is not evaluated in the simulations described in this work. 

B. Material State Mapping Algorithm 

Following crack growth and remeshing, state variables are mapped from the previous mesh to 
the current mesh using an inverse isoparametric mapping routine, as in [51]. Here, the scripts (n) 
and (n + 1) generically denote previous and current increments of crack growth, respectively. Lim et 
al. described the inverse isoparametric mapping technique for 2D elastic-plastic fracture simulations 
[52]. Implementation of the mapping algorithm consists of two high-level steps: (1) in the (n) mesh. 
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state variables stored at integration points are extrapolated to nodes using element shape functions 
and (2) displacements and state variables are transferred to either nodes or integration points in 
the (n + 1) mesh. The second step involves finding, for each point in the (n + 1) mesh, the natural 
coordinates (^, rf) of that point with respect to the element from the undeformed (n) mesh in which 
that point would spatially reside. The inverse problem becomes finding the natural coordinates 
(^, rf) that satisfy the known global coordinates: 

=SiVi(^,,?)Xy\ (5) 

where the subscript i ranges from one to the number of element nodes, are global nodal 

coordinates in the (n) mesh, are point or nodal coordinates in the (n + 1) mesh, and N are 

element shape functions evaluated at [C^r]). Once are known, nodal displacements and state 

variables, f/, can be transferred from the (n) mesh to the (n + 1) mesh in a forward manner, again 
using the element shape functions, N: 


7 ?) ( 6 ) 

Two levels of mapping are incorporated into the extended FRANC3D\NG and ABAQUS@ 
software framework. First, displacements are mapped onto the undeformed mesh following crack 
growth and remeshing. Second, a mapping function available in ABAQUS® is invoked to map the 
remaining state variables (e.g. stress, strain, plastic strain) onto the deformed mesh. When mapping 
material state between successive cracked configurations, it is critical that mesh refinement in regions 
of high gradients (e.g. near crack fronts) is sufficient to minimize solution diffusion, which occurs as 
a result of extrapolation, interpolation, and nodal averaging (if employed). This effect can become 
compounded as the crack growth simulation continues. After growing the crack and remeshing, 
the updated mesh model contains additional surface area due to crack extension, and equilibrium 
must be re-established before additional load is applied. During the equilibration procedure, global 
boundaries are held fixed and the new, traction-free crack surfaces are allowed to displace in response 
to surrounding fields, as shown in Fig. 10. 
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Fig. 10 Qualitative comparison of deformation and equivalent plastic strain field after mapping 
and subsequent equilibration. Images show face of 3D mesh. Deformation is not scaled. 


C. Validation Simulations 


Two stable tearing experiments are simulated to validate the EPFM framework for predicting 
crack growth and residual strength of relatively thin metallic structures. To illustrate the necessity 
of using an EPFM framework for predicting crack propagation and residual strength in such cases, 
the experiments are also simulated using an LEFM-based methodology. In the TEEM simulations, 
the material is modeled as linear-elastic, and crack growth is assumed to occur when an average value 
of Kj (and Ku for mixed-mode crack growth) along a crack front approximately equals fracture 
toughness of the material for any increment of crack length. 


1. Arcan Specimen 

A (modified) mixed-mode I/II Arcan fracture test [53] is simulated. Experimental details and 
results were provided in [54] and [55]. Drawings of the fracture specimen and modified load fixture 
are shown in Fig. 11. Curvilinear crack growth was induced in an A A 2024-T3 fracture specimen 
by applying monotonic load at a 30° angle relative to the fatigue precrack plane, as depicted in Fig. 
11 . 

The EE model of the load fixture and fracture specimen is constructed using ABAQUS@ and is 
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Fig. 11 Schematic of Arcan test fixture (left) and 2.3 mm thick fracture specimen (right). 
Load, P, can be applied at different pinholes in the fixture to induce mode I/II crack growth 
in the specimen. Drawings not to scale. All length dimensions in mm. 

depicted in Fig. 12. The FE model contains an initial crack of length 6.35 mm and does not simulate 
the fatigue precracking process. A mesh refinement study is first conducted by simulating a pure 
mode I Arcan test (0° loading angle). Results from the study reveal that, for the initial crack, the 
applied load at CTDcrit varies by less than 0.25% when the element length nearest the crack front is 
0.33 mm or smaller. This converged level of mesh refinement is used for the mixed-mode Arcan model 
presented here (30° loading angle). The specimen and fixture are assumed to be perfectly bonded 
such that coincident nodes are merged at the specimen-fixture interface. The fixture is modeled using 
6000 C3D10 elements and remains unchanged throughout the simulation. The fracture specimen 
sub-region is subject to geometry and mesh updating within FRANC3D\NG. Depending on the 
crack length, the specimen comprises between 9000 and 25,000 quadratic elements, which include 
a standard rosette of C3D15, C3D20, and pyramid (collapsed C3D20) elements surrounding the 
crack front (see [42] for details). The bulk of the sub-region mesh comprises C3D10 elements. 

Material properties for the 15-5PH stainless steel fixture are E = 207 GPa and v = 0.3. The 
fixture is assumed to remain elastic during loading. Material properties for the AA 2024-T3 fracture 
specimen are E = 71.4 GPa, u = 0.33, and ay = 345 MPa. The strain hardening curve for the 
specimen is provided in Fig. 13. A von Mises yield criterion with isotropic hardening is assumed. 
For the LEFM simulation, the specimen is modeled as linear-elastic with Kjc = 37 MPa^m for 
A A 2024-T3 in the LT orientation [56]. 

Crack growth occurs in the LEFM simulation when the average (AT/, Kn) along the crack front 
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Fig. 12 3D FE model of modified Arcan test set-up, including load fixture and fracture spec- 
imen. Line traction, P, is applied at 30° from mode I axis (^/-axis). 



Fig. 13 Strain hardening curve determined from uniaxial tension tests for AA 2024-T3 in LT 
orientation. Courtesy NASA Langley Research Center. Similar curves were used in [22] and 

[23]. 


reaches a critical combination. The critical combination is determined using the mixed-mode failure 
envelope described in II B, where Kuc is taken to be 10% less than Kjc. 

In the EPFM simulation, important observations from a number of studies inform the selection 
of CTDcrit and d used to predict crack growth. First, from [6, 54, 55], for 2.3 mm thick AA 2024-T3 
specimens precracked in the L-T orientation, an average constant critical value of CTDcrit= 0-1 nim 
was observed, where CTD was measured on the specimen face at d = 1 mm. This critical value was 
observed for a range of specimens exhibiting mode I dominant crack growth, which includes the 30° 
Arcan test [54, 55]. Scatter among the measurements was typically ±0.02 mm. Second, significant 
crack front tunneling was observed in the fracture specimens, especially during initial crack extension 
[6]. Third, a recent study by Lan et al. suggested that modeling tunneled cracks using straight crack 
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fronts can lead to over-estimation of load versus crack extension predictions [50]. This is because 
higher levels of constraint along a crack front, which can induce crack tunneling, effectively decrease 
fracture toughness. Lan et al. noted that surface-measured CTD values can be 24% larger for cases 
with crack tunneling than cases without. In the Arcan simulation, crack front tunneling and slant 
crack growth are not modeled. Thus, using the surface-measured value of CTDcrit= 0.1 mm from 
[6, 54, 55] as the crack growth criterion in the Arcan simulation, an over-prediction of load versus 
crack extension is expected. Assuming that the over-prediction is proportional to the difference in 
surface-measured CTDcrit values between straight and tunneled crack fronts, CTDcrit is taken to be 
0.08 mm, which is ^ 24% less than the observed average surface-measured value of CTDcrit= 0.1 mm 
and also corresponds to the lower bound of experimental scatter. This adjusted critical value used 
in simulation should account for some expected over-prediction of residual strength. 

Simulated crack growth occurs in increments of approximately 1 mm, which satisfies convergence 
results from similar simulations [22]. Propagation direction is predicted according to the maximum 
tangential stress theory [16], though a 2D CTD-based directional criterion [19] could also be used. 

Simulation proceeds by applying displacement in the direction indicated in Fig. 12, which cor- 
responds to a 30° loading angle. In the EPFM simulation, after each increment of, inclusively, 
crack growth, remeshing, and material state mapping (see subsection IIIB), all nodes with applied 
boundary conditions are held fixed while the model is brought into equilibrium. Then, additional 
displacement is applied, and the simulation continues until maximum load, or residual strength, is 
attained. 

Applied load versus crack extension is plotted in Fig. 14 from the Arcan experiment and sim- 
ulation. Both LEFM and EPFM simulation results are plotted. Compared to experiment, residual 
strength is 6.8% higher using the EPFM framework and 10.8% lower using the LEFM method. 
In the LEFM simulation, maximum load occurs at the first increment of crack growth; whereas, 
maximum load occurs in the EPFM simulation after a small amount of crack growth, consistent 
with experiment. The load to initiate crack extension is approximately the same in both EPFM and 
LEFM simulations since there are no residual stresses in the model at da=0. However, unlike the 
LEFM simulation, the EPFM simulation is able to capture the shape of the load versus crack exten- 


27 



sion curve. In other words, including plasticity effects in the model requires an increase in applied 
load, initially, to drive crack extension before reaching the residual strength of the specimen. 

Although the EPFM framework captures the shape of the load versus crack extension curve, it 
slightly over-predicts the maximum load. The over-prediction is likely related to modeling the crack 
front, which tunnels or thumbnails in reality, as remaining straight throughout the simulation. Given 
a constraint-dependent CTDcritj curvature along the crack front could be predicted using a point- 
by-point evaluation technique, in which case it is suspected that the load versus crack extension 
curve would be predicted even more accurately. Another potential contribution to the discrepancy 
is the numerical tolerance for satisfying CTDcritj which in this case is set at 2%. 
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Fig. 14 Applied load versus crack extension, da, from the 30° Arcan experiment and sim- 
ulations. Insets show snapshots of deformed mesh at various crack increments. Complete 
simulation can be viewed at www.cfg.cornell.edu 


Figure 15 shows the simulated and experimental curvilinear crack trajectory as viewed from 
the specimen free surface. There is a slight deviation in the actual crack trajectory during the 
first 6.35 mm of crack length due to fatigue precracking processes, which are not modeled in the 
simulation. Rather, the fatigue precrack is simply modeled as a perfectly planar initial crack in the 
simulation. The deviation appears to have little effect on the predicted crack path thereafter. 
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Fig. 15 Comparison of experimental and simulated curvilinear crack paths in the 30° Arcan 
fracture test. Inset shows reference coordinate system. 

2. Integrally- stiffened Panel 

A stable tearing test of an ISP machined from a lower wing-skin aluminum alloy, C433-T39, is 
also simulated. The test was conducted at Alcoa Technical Center. Test details, data, and results 
were overviewed in [57] and have also been provided to the authors by Alcoa Technical Center. 
Relevant details are described next, and additional details from the test program can be found in 
the Appendix. 

Dimensions of the panel are shown in Fig. 16(a). An initial two-bay saw cut of length ~ 2.54 cm 
was made at mid-height to completely sever the middle stiffener. The initial cut was then propagated 
under fatigue loading until both crack fronts were 2.54 cm short of reaching the intact stiffeners 
(2a ~ 24.1 cm). The panel was then loaded monotonically in uniaxial tension until failure occurred 
by unstable crack growth. Crack front branching was observed, where an initial crack propagating 
toward an intact stiffener eventually split into two distinct cracks, one continuing into the adjacent 
bay and one propagating in the 2 ; direction within the stiffener. Photographs of the test panel with 
views of crack branching are provided in Fig. 17. 

A 3D FE model of the entire panel is constructed using ABAQUS@ [1]. The FE model contains 
an initial crack of total length 24.1 cm, which corresponds to the fatigue crack length just prior to 
conducting the residual strength test. The FE model with initial crack is shown in Fig. 16(b). The 
mesh region that remains unchanged throughout the tearing simulation is modeled using 56 C3D15 
and 9,400 C3D20R elements. A 38.5 x 12.7 cm^ sub-region centered in the panel is subject to 
geometry and mesh updating within FRANC3D\NG. Depending on crack length, the sub-region 
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comprises between 27,000 and 95,000 quadratic elements, including a bulk of (731)10 elements and 


a standard rosette of (73D15, (73D20, and pyramid (collapsed (73D20) elements surrounding the 
crack front (see [42] for details). The mesh interface between the sub- and global regions is coherent, 
obviating the enforcement of a constraint. 


3.8 



10.2 15.2 
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Fig. 16 (a) Schematic (not to scale) from [57] showing dimensions of symmetric ISP tested at 
Alcoa Technical Center. Isoparametric view of full panel indicates final fatigue crack length 
2a^, and cross-section view in plane of the crack shows where stiffener is completely severed. 
All dimensions are in cm. (b) Corresponding 3D FE model of ISP. Initial crack severs middle 
stiffener. Traction, P, is applied uniaxially in the y direction. 


Mesh refinement near each crack front is dictated, to some extent, by crack front location with 
respect to stiffener. As the crack initially propagates through the skin, element lengths nearest 
both crack fronts are 0.5 mm. When the crack fronts are within 1.5 mm of the 90° skin-stiffener 
junctions, however, near-crack front elements must decrease in size. This is because the rosette 


30 



Outboard view 

Fig. 17 Photographs of ISP with central two-bay crack from the Alcoa test program [57]. Full 
panel in load frame (left) and angled views of crack front branching into stiffener (top right) 
then exiting the stiffener (bottom right). 


template of elements surrounding either crack front consists of three rings of equi-length elements. 
In order to accommodate (i.e. facilitate remeshing) a full rosette of elements within the proximity 
of the discontinuous geometry, the size of the template elements must decrease. 

The thickened grip ends of the panel are modeled as linear-elastic with an elastic modulus 
approximately five times greater than that of C433-T39. The rest of the panel is assigned C433-T39 
material properties: E = 71 A GPa, u = 0.3, and ay = 455 MPa [58]. The strain hardening curve 
used for C433-T39 is provided in Fig. 18. A von Mises yield criterion with isotropic hardening is 
assumed. For the LEFM simulation, the panel is modeled as linear-elastic with Kjc = 50 MF^i^/rn 
for C433-T39 [58]. 



Fig. 18 Strain hardening curve determined from uniaxial tension tests for C433-T39 in LT 
orientation [57]. 


Boundary conditions are applied to simulate actual loading in the panel. Nodes on the bottom 
face of the lower grip end are fixed in the y direction. Displacement is applied in the y direction at 
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nodes on the top face of the upper grip end. Additionally, nodes along the same top and bottom 
grip end faces are fixed in the x and 2 ; directions. For the EPFM simulation, after each increment 
of, inclusively, crack growth, remeshing, and material state mapping (see subsection IIIB), all 
nodes with applied boundary conditions are held fixed while the model is brought into equilibrium 
before applying additional displacement. Additionally, based on preliminary simulation results, the 
entire back (zmin) face is artificially fixed after mapping and during the equilibration procedure. 
This is because resonance in the 2 ; direction is observed with increased crack growth otherwise. 
The resonance occurs when mapped tensile and compressive stresses in the faces of the panel are 
artificially reversed during equilibration of the mapped solution. The additional boundary condition 
is, however, removed after the equilibration procedure so that 2 ; displacement is allowed during the 
subsequent loading step. 

Crack growth occurs in the LEFM simulation when the average Kj value along either crack 
front reaches Kjc. A mixed-mode failure criterion is unnecessary, as Ku and Km are negligible 
(i.e. <2.5% of Ki for all crack growth increments). 

For the EPFM simulation, CTDcrit was calibrated at NASA Langley Research Center from a 
middle-crack tension (MT) test of the same material (C433-T39) and thickness as the ISP [59]. In 
that work, 3D FE simulations of the MT test revealed that simulated load versus crack extension 
matched experimental data when the mode I opening angle midway along the crack front at d = 
1.02 mm reached a critical value of 6.5°. This angle corresponds to CTDcrit through the relation 

tan(6.5°) = (7) 

The same criterion is applied in the EPFM simulation by specifying for both crack fronts that 
CTDcrit must attain a value of 0.116 mm at d = 1.02 mm behind the crack front and that the 
criterion be evaluated midway along either crack front. 

Crack growth occurs in increments of 1.15 mm (about 15% of the skin thickness), which is 
selected to be approximately the same as that implemented by Seshadri et al. in a similar simulation 
[59]. Straight crack fronts are enforced during crack growth in the skin of ISP. Upon entering the 
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intact stiffener, a crack front is evolved such that (1) a realistic, arbitrary crack front profile is 
represented (though the actual evolving crack front profile was not monitored during experiment) 
and (2) the new crack front profile has relatively smooth curvature to facilitate remeshing. A cross- 
section view of the panel in Fig. 19 shows different stages of simulated crack front evolution, from 
lead crack growth in the skin, to transition crack growth within the stiffener, to complete branching. 
The simulation proceeds as depicted in Fig. 8 until both initial crack fronts completely branch and 
Pmax is attained. 




e) _ f) 

Fig. 19 Cross-sectional views of ISP mesh model taken at the crack plane and magnified at 
one stiffener. A thickened red line is overlaid along the crack front (s) at each step of crack 
growth. Views show lead bay crack before entering stiffener (a); transition crack evolution 
within stiffener (b,c,d); and complete branching into two distinct crack fronts (e,f). 


Evaluation of the CTD criterion becomes nontrivial as a crack front transitions within the stiff- 
ener (i.e. while a crack front is within the stiffener but has not completely branched). Constraint 
effects introduced by the stiffener on the unsymmetric crack front profile, along with slight 2 ; dis- 
placement near the cracked region, lead to nonuniform and unsymmetric CTD values along the 
crack front. Evaluating the CTD criterion at only one point along the crack front becomes ambigu- 
ous to implement numerically and less representative, physically, of 3D crack growth behavior. A 
simple and efficient approach to address these issues is to compare CTDcrit to an average of CTD 
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values along the crack front. For transition crack growth in the stiffener, the middle third section of 
CTD points along the crack front are averaged and evaluated to predict crack propagation. Once 
the crack fully branches, the CTD criterion is again evaluated midway along each crack front. 

Predicting crack front evolution for transition crack growth within the stiffener is not currently 
a fully-automated process. In this work, the straight portion of the crack front (see Fig. 19) is 
predicted using FRANC3D\NG and the curved portion of the front is specified by manually adding 
crack front points according to the two considerations described in the previous paragraph. Given a 
constraint-dependent CTDcrit relationship, a point-wise CTD criterion could be evaluated to evolve 
the arbitrarily-shaped crack front automatically. 

The consequence of using an LEFM versus EPFM simulation to determine residual strength is 
quite obvious in the ISP case. Figure 20 shows load versus crack extension for both LEFM and EPFM 
simulations. Experimental load versus crack extension was not recorded during the tests; however, 
maximum applied load is plotted for two ISP tests of the same material and loading conditions. As 
in the Arcan simulation, the load to initiate crack extension is similar using either EPFM or LEFM 
methods since there are no residual stresses in the model at da=0. Following initiation, however, the 
EPFM simulation predicts that the applied load must be increased to maintain crack propagation. 
The necessary increase in applied load occurs since a significant amount of energy in the system is 
dissipated through plastic deformation. This effect cannot be predicted by the LEFM simulation 
since plasticity effects are not modeled. As a result, much of the energy in the ISP for the LEFM 
simulation must be dissipated through creation of new fracture surface area, which means less load 
is required to drive crack growth in the LEFM simulation than in the EPFM simulation. Using the 
EPFM framework, residual strength is determined within 2% of experimental average of the two 
tests. On the other hand, the LEFM method underpredicts the average residual strength by 64%. 
Although the LEFM simulation does predict an increase in applied load at the stiffener junction 
due to geometrical effects, the increase is negligible compared to that due to plastic deformation. 

From the EPFM simulation, equivalent plastic strain field evolution in the ISP is depicted 
in Fig. 21 for the first and final crack steps. Accumulation of plastic strain in the wakes of the 
advancing crack fronts is relatively significant, extending from initial to final crack front locations. 
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Fig. 20 Applied load (traction P, Fig. 16(b), integrated over applied area) versus half-crack 
extension, da, from ISP simulation. Maximum applied load is indicated for two corresponding 
tests conducted at Alcoa Technical Center. Shaded region indicates initially intact stiffener. 


The general shape of 45° contour lobes at da=0 mm is maintained at da =44 mm both for the lead 
crack front extending into adjacent bay and for the crack front propagating in the 2 ; direction within 
the stiffener. At da=44 mm, the equivalent plastic strain in each stiffener extends in the direction 
of both contour lobes to the stiffener boundary. The consistent contour lobe shapes indicate that, 
despite increased 2 ; displacement as the crack propagates and severs initially intact stiffeners, both 
lead and branched crack fronts remain locally mode I dominant throughout tearing. 

Finally, as evident in Fig. 21 for da=44 mm, the mapping procedure inevitably leads to imper- 
fections in the fields due to diffusion of the FE solution, which occurs in regions of high gradients 
from repeated extrapolation and interpolation procedures, see IIIB. If mapping errors significantly 
affect crack growth predictions, mesh refinement should mitigate this effect. 


IV. Conclusions 

A surrogate model methodology and 3D elastic-plastic fracture simulation toolset have been 
presented, which enable accurate residual strength prediction of damaged structures in real time. 
The methodology and toolset are particularly useful for scenarios involving metallic aircraft struc- 
tures subject to discrete-source damage during flight. An accurate prediction of structural residual 
strength in these scenarios could aid in avoidance of catastrophic crack growth and subsequent 
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Fig. 21 Magnified views of simulated crack growth in ISP at half-crack extensions da=0 mm 
(top) and da=44 mm (bottom). Contours show evolution of equivalent plastic strain fields 
with crack growth. Deformation is not scaled. FE mesh is not shown for better contour 
visualization. Complete simulation can be viewed at www.cfg.cornell.edu. 


structural failure. 

The surrogate model methodology relies on offline numerical fracture simulations to obtain a 
set of data points describing residual strength as a function of discrete-source damage parameters. 
Strictly for illustration, a NN has been constructed as a surrogate model for predicting residual 
strength of a representative wing sub-structure subject to discrete-source damage. In the illustra- 
tion, offline residual strength values have been determined using computationally efflcient LEFM 
approximations. Subsequently, the consequences of using LEFM approximations for determining 
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residual strength of damaged metallic structures have been shown, and an EPFM framework to 
accurately determine residual strength using high-fidelity, 3D, elastic-plastic tearing simulations has 
been described. For an aluminum-alloy, integrally-stiffened panel exhibiting crack branching, resid- 
ual strength is predicted within 2% of experiment using an EPFM simulation and is underpredicted 
by 64% using an LEFM simulation. 

The more general and rigorous elastic-plastic tearing framework should be used to generate ac- 
curate residual strength training data, especially for cases involving discrete-source damage. Also, 
the FE model for the structure of interest should include enough detail to fully capture the rela- 
tionship between a particular global loading state and onset of unstable crack growth. Furthermore, 
damage should be parameterized by taking into account onboard sensor characterization capability 
and resolution. With these considerations in mind, the general surrogate model methodology cou- 
pled with the EPFM simulation framework presented in this work provides a means of achieving 
more resilient and adaptive aircraft control. 
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VII. Appendix 

A. Python@ Script to Evaluate CTD Criterion During Nonlinear Finite Element Analysis 
Using ABAQUS® 

The script CTDjohControlTerminal.py executes an FE analysis using ABAQUS® and period- 
ically interrupts the analysis to compute CTD values along one or multiple crack fronts and to 
subsequently evaluate a user-specified CTDcrit criterion. If CTDcrit is satisfied within the given 
tolerance, the script will terminate the FE analysis and echo to the command window the load 
increment when the criterion is satisfied. If there are multiple cracks in the model, the critical 
crack front will be identified and echoed to the command window. The script will also print the 
file CTDResults.txt to the current working directory. The file includes CTD values (decoupled into 
three modes) at each point along each crack front for every FE analysis interruption. The script 
includes function calls to the object libraries MeshTools and Vec3D, which are not included here. 
However, the functionality of each routine should be somewhat clear from the name of the function 
and comments in the script. 

For the first step of crack growth, assuming the initial state is undeformed, CTD job Control, py 
should be located in a current working directory and should be executed from the MS-DOS command 
window using a command like: “...>abaqus python CTDJobControl.py”. The user will then be 
prompted for a series of input. Successful execution requires that the following files be located in the 
current working directory along with the script: (1) *.fdb (generated by FRANC3D\NG following 
crack insertion and remeshing) and (2) the deformed *.inp (only the model to be submitted for FE 
analysis). Additionally, the undeformed *.inp file (the undeformed local mesh model in a global- 
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local analysis) should be included in a sub-directory called Undeformed within the current working 
directory. 

For subsequent steps of crack growth, a similar script, CTDjohControlCallpy, is called auto- 
matically following deformation mapping. Unlike CTDjohControl.py, CTDjohControlCall.py does 
not require the user to relocate files or execute a script to start the next FE analysis. While the 
two scripts differ slightly in how they are executed, the main functions are the same. Only CTDjoh- 
ControlCall.py is included here, and electronic versions of both scripts are available for download 
at www.efg.eornell.edu. 


? ? ? 

CTD j obControlCall . py 

Python script to control ABAQUS nonlinear FEA while 
monitoring CTD at points located distance d behind crack, 
written by ADS (October, 2009) 
modified March, 2010 

***if executing manually, use: "...>abaqus python CTDj obControlCall .py -- 
def InpFilePath undef InpFilePath fdbFilePath oldJobodbPath"*** 

Otherwsie the script will be executed automatically after mapping 
deformation using runMapScript .py 

? ? ? 

import os, sys, time 

from sys import argv, exit 

import odbAccess as oa 

from abaqusConstants import * 

import VecSD, MeshTools 

import math 

from math import sqrt 

# ************************** define helper functions ******************************* 

def askForFloat (number) : 

# Request float from user 
response = float (raw_input (number) ) 
if response < 0 : 

print "Requested value must be positive. Exiting." 
exit (0) 
else : 

return response 
def askForInt (number) : 
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# Request int from user 
try: 

response = int (raw_input (number) ) 
except : 

print "Response must be a positive integer!" 
exit (0) 

if response < 0 : 

print "Requested values must be a positive integer. Exiting." 
exit (0) 
else : 

return response 

def askForPath (prompt , extension) : 

# Request path from user 
response = raw_ input (prompt) 
while not os .path. exists (response) : 

print "Specified path could not be located." 
print "Check that path exists and try again." 
response = raw_ input (prompt) 
while not response . split (’.’) [-1] == extension: 

print "\n**File path should have extension ’Xs’**" % extension 
response = raw_ input (prompt) 
return response 

def changeDirectory (pathToFile) : 
path = str() 

tmp = pathToFile . split ("/") 
for i in range (len (tmp) -1) : 
path+=tmp [i] 
path+="/" 

print "\nThe cwd has been changed to: \n\t7oS" % path 
os . chdir (path) 

def getFdb(path) : 
fdb = str() 

# First try splitting on "." 
tmp = path. split (". ") 
tmp[-l] = ’.fdb’ 

for i in range (len (tmp) ) : 
f db+=tmp [i] 

# If doesn’t exist, try splitting on "_" 
if not os .path. exists(fdb) : 

fdb = str() 
tmp = path. split ("_") 
tmp[-l] = ’.fdb’ 
for i in range (len (tmp) ) : 
f db+=tmp [i] 

if not os .path. exists (fdb) : 

print "The file Xs cannot be located." % fdb 
exit (0) 

print "\nGetting Xs" % fdb 
return fdb 
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def askYesNo (question) : 

# Request yes/no response from user 

# (Returns "Y" or "N") 
response = raw_ input (question) 

while response .upper 0 [0] not in ("Y","N"): 
response = raw_ input (question) 
if response .upper 0 [0] not in ("Y","N"): 
print "Wrong answer! Try again" 
return response .upper () [0] 

def isint (value) : 

# Check if value is an integer 
try: 

int (value) 
return True 
except : 

return False 

def getGroupList (meshObj , listName , type) : 

list = meshObj .GetGroupInfo (listName, type) 
if len(list) == 0: 

print "***No members found in group: Xs ***" % listName 
return list 

def sortDictValsByKey (diet) : 

# Returns a list of values sorted by keys 
sorted_list = [] 

keys = diet. keys 0 
keys . sort 0 

for i in keys: 

sorted_list . append (diet [i] ) 

return sorted_list 

def printTimeStamp(start_time) : 

# Print time stamp 
end_time = time.timeO 

mins = int ((end_time-start_time)/60) 
secs = round( (end_time-start_time)7o60,3) 
resultsFile = ’CTDResults.txt’ 
file = open(resultsFile , "a") 

file . write ( ’ \n* **************************************** ********** ’ ) 
file, write ( ’\nTotal time for run: 7oi:7og’ 7o (mins, secs)) 
file . write ( ’ \n* **************************************** ********** ’ ) 
file. closeO 

# **************** define job/analysis interaction functions *********************** 

def runAnalysis (d , CTDc , eval , oldJob ,mapSoln ,undef InpF ile , inpFile , f dbFile) : 
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# ************************* runAnalyis SUB-FUNCTIONS *************************** 
def submit Job 0 : 

# Submits ABAQUS job and suspends it after data written . sta file 
print "\nSubmitting ABAQUS job ’7oS’..." % jobName 

if oldJob in "None": 

os . system (’ abaqus job=’+jobName) 
else : 

print "Old job specified" 

os . system ( ’abaqus job=’+jobName+’ oldjob=’+oldJob) 
while not os .path. isf ile(staFile) : # Wait for .sta file to be created 

continue 

# A hack at monitoring the sta file... 
fileSize = os .path.getsize(staFile) 
print fileSize 

f ileSize_update = fileSize 
while fileSize == f ileSize_update : 

f ileSize_update = os .path.getsize(staFile) 
else : 

print f ileSize_update 

os . system (’ abaqus job=’+jobName+’ suspend’) 

def resumeJobO : 

# Resume ABAQUS job 

os . system ( ’abaqus job=’+jobName+’ resume’) 

# Monitor the sta file 

fileSize = os .path.getsize(staFile) 
print fileSize 
f ileSize_update = fileSize 
while fileSize == f ileSize_update : 

f ileSize_update = os .path.getsize(staFile) 
else : 

print f ileSize_update 

os . system ( ’abaqus job=’+jobName+’ suspend’) 

def openOdbC jobName) : 

odbFile = j obName+ ’ . odb ’ 

# Try opening the odb file 
ct = 1 

while ct < 11: 
try: 

odb=oa . openOdb (path=odbFile , readOnly=TRUE) 
print "Opening the odb file..." 
return odb 
except : 

print "ERROR: Unable to open the specified odb Xs ! " % odbFile 
ct+=l 
exit (0) 

def getCurrentStatus 0 : 

fd = open(staFile , "r") 
lines = fd. readlines 0 
i = 1 

if len(lines [-i] . split 0 ) == 0: 
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i+=l 

while isint (lines [-i] . split 0 [0] ) == False : 
i+=l 

while len(lines [-i] . split 0 ) == 0: 
i+=l 

curr_step = int (lines [-i] . split () [0] ) 
curr_inc = int (lines [-i] . split () [1] ) 


return curr_step, curr_inc 

def previous Increment (f rac_status_per_f ront , front , eval , curr_inc , curr_step , \ 

interest ing_nodes ,mapped_disp ,pts_per_f ront , \ 
main_side_elems ,mate_side_elems , odbStep) : 

# Evaluate CTD values at successive previous increments 

# until no values are supercritical, 
while frac_status_per_f ront [front] == 2: 

if curr_inc == 1: 

print "\nCTD(A) exceeds critical at crack front Xi on load \ 
increment 1" % front 

print "Decrease initial inc size! See CTD results file for details" 
os . system ( ’abaqus job=’+jobName+’ terminate’) 
exit (0) 
else : 

print "\nStepping back one increment..." 
curr_inc-=l 

relative_disp = disps_per_inc [curr_inc] 
if not mapped_disp: 

abs_disp = relative_disp 
else : 

abs_disp = getAbsDisps(mapped_disp,relative_disp) 
updateDisplacements (mesh, abs_disp) 

f lags_per_front = computeAndWriteCTDResults(mesh,tol,d,CTDc,\ 
inpF ile , f dbF ile , j obName , odbStep , curr_inc , \ 
pts_per_f ront ,main_side_elems ,mate_side_elems) 
frac_status_per_f ront = evaluateFractureCriterion(f lags_per_front ,eval) 
print "\nCrack front Xi was supercritical at\n7oS increment 7oi" \ 

7o (front , odbStep, curr_inc+l) 

print "May need to decrease FEA increment! See results file." 


# sie***************************** runAnalysis 

# Initializing required inputs... 

jobName = inpFile . split ("/") [-1] . split (".") [0] 
staFile = jobName+" . sta" # 

base = undef InpFile . split (".") [0] # 

mesh = MeshTools .MeshTools(base, "INP") 
inc_interval = 1 

i = inc_interval # 

disps_per_inc = {} # 




# Job name corresponds to . inp file 
Status file corresponds to job 
Undeformed .inp file in \Undeformed subdir. 


_per_in 
disps_per_inc [0] = 0 
f lags_per_front = {} 
eq_flag = 1 


Access the odb at every ith FEA increment 
A diet containing incremental displacements 
Initialize the dictionary 


# Flag to indicate equilibration step 


# Calling CTD initializing functions... 

# Get mapped displacements from mapped_disps.txt (if any) 
mapped_disp = get InitialDisps (mapSoln, inpFile) 
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# Get main/mate side node groups 

(main_side_nodes , mate_side_nodes) = getCrackFaceNodes (mesh) 

# Main/mate side elements along crack faces 

(main_side_elems , mate_side_elems) = getCrackFaceElems (mesh,main_side_nodes , \ 

mate_side_nodes) 

# Nodes belonging to main/mate side elements (only need to update disps 

# for these nodes during the analysis for computing CTD values) 

interest ing_nodes = getCrackFaceElemNodes (mesh,main_side_elems ,mate_side_elems) 

# Determine points behind crack front (s) where CTDs computed ("CTD points") 
pts_per_f ront = getPoints (main_side_elems ,mate_side_elems ,main_side_nodes , \ 

fdbFile ,mesh,d) 

# Set tolerance distance 

tol = setTol(mesh,main_side_nodes) 

# JOB CONTROL: 

# Initial job submission, suspension, and CTD computation 
start_time = time. time () 

submit Job 0 

(curr_step, curr_inc) = getCurrentStatus () 

# Keep track of current step for disps_per_inc dictionary 

# (clear the dictionary when we start a new step) 
step = curr_step 

# In some cases, increment info may not be written to odb yet, 

# in which case disps_per_inc (and f lags_per_front) remains empty, 
while not f lags_per_front : 

print "Waiting for displacements to be written to odb..." 

odb = openOdb(jobName) # The odb must be closed then re-opened if 

# we need information that has not been written 

# to the odb at the current time of access 
(curr_step, curr_inc) = getCurrentStatus () 

(disps_per_inc , odbStep) = getRelDispsAtInterestingNodes (disps_per_inc , \ 

interesting_nodes , j obName , curr_step , curr_inc , odb) 

# Raise flag if done with equilibration step 
if odbStep .upper 0 [0 : 3] != "EQU" : 

eq_flag = 0 


try: 

relative_disp = disps_per_inc [curr_inc] # Disps available at current inc? 
if not mapped_disp: # If no initial mapped disps, 

# then relative_disp=abs_disp 
abs_disp = relative_disp 
else : 

abs_disp = getAbsDisps (mapped_disp ,relative_disp) 
updateDisplacements (mesh , abs_disp) 

f lags_per_f ront = computeAndWriteCTDResults (mesh, tol ,d,CTDc , \ 

inpF ile , fdbFile , j obName , \ 
odbStep, curr_inc ,pts_per_f ront , \ 
main_side_elems ,mate_side_elems) 

except : 

odb . close 0 
resumeJobO 
continue 
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# While the script is running, perform actions based on CTD flags returned. 

# frac_status_per_front --> A diet where key is crack front id and value is: 

# 0 --> ’ subcritical ’ , continue 

# 1 --> critical, extend crack 

# 2 --> ’supercritical’, back-up or decrease FEA inc 
running = True 

while running: 

if eq_flag == 0: 

# Evaluate the fracture criterion, sending to the function the flags 

# for all CTD (A) points at each crack front and the evaluation criterion 

# to be used 

frac_status_per_front = evaluateFractureCriterion(f lags_per_f ront , eval) 

# First loop through all crack fronts and check for any that are 

# supercritical : (if there are, then the analysis will terminate as soon 

# as this is not the case or inform user to decrease load increment) 
for front in f rac_status_per_front : 

if frac_status_per_f ront [front] == 2: 

previousincrement (f rac_status_per_f ront , front , eval , curr_inc , \ 
curr_step , interesting_nodes ,mapped_disp , \ 
pts_per_f ront ,main_side_elems ,mate_side_elems , \ 
odbStep) 

os . system ( ’abaqus job=’+jobName+’ terminate’) 

printTimeStamp(start_time) 

exit (0) 

# Next loop through all crack fronts and check for any that are critical: 

# (if there are, then analysis will terminate at current step/increment) 
critical = 0 

for front in f rac_status_per_front : 

if frac_status_per_f ront [front] == 1: 

print "\n***Crack front Xi is within tolerance of critical***"\ 

7o front 
critical = 1 
if critical == 1: 

print "\nExtend critical crack (s) at increment Xd of step 7oS"\ 

7o (curr_inc , odbStep) 

os . system ( ’abaqus job=’+jobName+’ terminate’) 

printTimeStamp(start_time) 

exit (0) 

# Finally, if analysis has not terminated, all crack fronts are subcritical 
odb. close 0 # Close odb to update add’l analysis increments 

resumeJobO 

(curr_step, curr_inc) = getCurrentStatus () 

# If we started a new step, clear disps_per_inc dictionary and start over 
if curr_step > step: 

disps_per_inc . clear () 
disps_per_inc [0] = 0 
step = curr_step 
odb = openOdb(jobName) 

(disps_per_inc , odbStep) = getRelDispsAtInterestingNodes (disps_per_inc , \ 

interesting_nodes , j obName , \ 
curr_step , curr_inc , odb) 

# Raise flag if done with equilibration step 
if odbStep .upper 0 [0 : 3] != "EQU" : 
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eq_flag = 0 


relative_disp = disps_per_inc [curr_inc] 
if not mapped_disp : 

abs_disp = relative_disp 
else : 

abs_disp = getAbsDisps (mapped_disp ,relative_disp) 
updateDisplacements (mesh,abs_disp) 

f lags_per_f ront = computeAndWr iteCTDResult s (mesh , tol , d , CTDc , inpFile , \ 

f dbF ile , j obName , odbStep , curr_inc , \ 
pts_per_f ront ,main_side_elems , \ 
mate_side_elems) 


# * * * * * * * * * * * * * * * * * * * * * * * * * * DEFINE CTD FUNCTIONS ********************************* 
def getInitialDisps (mapSoln,f ilePath) : 

# Retrieve mapped displacements from mapped_disp.txt. Return empty diet if none. 
mapped_disp = {} # Diet of mapped (initial) nodal displacements 

if mapSoln == "N" : 

print "***Returning empty mapped_disp diet***" 
return mapped_disp 

elif mapSoln == "Y" : 

changeDirectory (f ilePath) 
disp_file = ’mapped_disp.txt’ 

# If the disp.txt file doesn’t exist, exit 
if not os .path. exists (disp_f ile) : 

print "***Could not locate mapped_disp.txt from previous analysis!***" 
print "(Make sure the file is located with the previous mesh file)" 
exit (0) 

file = open("mapped_disp. txt" , "r") 
buff = file .readline 0 
while buff : 

data = buff. split 0 
if data [0] .upper 0 [0:4] == "DISP": 
buff = file. readline 0 
data = buff. split 0 
while buff : 

nid = int(data[0]) 
tmp = [] 

for i in range (3): 

tmp. append (float (data[i+l] ) ) 
mapped_disp [nid] = tmp 
buff = file .readline 0 
data = buff. split 0 
buff = file .readline 0 

file . close 0 

print "***Retrieved initial displacements for Xd nodes" % len(mapped_disp) 
return mapped_disp 
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def getCrackFaceNodes (mesh) : 

main_side_nodes = [] 
mate_side_nodes = [] 

try: 

main_side_nodes = getGroupList (mesh, "main_side_nodes" , "node") 
mate_side_nodes = getGroupList (mesh, "mate_side_nodes" , "node") 
except : 
try: 

main_side_nodes = getGroupList (mesh, "all_main_side_nodes" , "node") 
mate_side_nodes = getGroupList (mesh, "all_mate_side_nodes" , "node") 
except : 

print "***No main/mate side nodes found for previous mesh!***" 
exit (0) 

return main_side_nodes , mate_side_nodes 

def getCrackFaceElems (mesh,main_side_nodes ,mate_side_nodes) : 

# Make sets of main/mate side crack face elements 

print "\nCollecting elements along the main and mate side crack face..." 
main_side_elems = set() 
mate_side_elems = set() 

for nid in main_side_nodes : 

elist = mesh. Get AdjacentElems (nid) 

for eid in elist: main_side_elems . add(eid) 

for nid in mate_side_nodes : 

elist = mesh. Get AdjacentElems (nid) 

for eid in elist: mate_side_elems . add(eid) 

return main_side_elems , mate_side_elems 

def getCrackFaceElemNodes (mesh,main_side_elems ,mate_side_elems) : 

# Gather only nodes of interest for updating displacements 
nodes = [] 

for eid in main_side_elems : 

nids = mesh. GetElemInfo (eid) 
for n in nids: 

if n not in nodes : 
nodes . append (n) 

for eid in mate_side_elems : 

nids = mesh. GetElemInfo (eid) 
for n in nids: 

if n not in nodes : 
nodes . append (n) 

return nodes 

def getPoints (main_side_elems ,mate_side_elems ,main_side_nodes , f dbFile ,mesh , d) 
print "\nDetermining coordinates of CTD(A) points Xg behind \ 
crack front ..." % d 
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# First, create a dictionary of crack fronts as keys and 

# lists of crack front node id’s as values (from the fdb) 
front_nodes = {} 

nodes = [] 
f=0 


fdb = open(f dbFile , ’r ’ ) 

buff = fdb .readline 0 
while buff : 

vals = buff. split 0 
if vals [0] .upper 0 == "NUM.FRONTS : " : 
nfronts = int(vals[l]) 
break 

buff = fdb.readlineO 

buff = fdb.readlineO 
while len(front_nodes) < nfronts: 
buff = fdb.readlineO 
vals = buff. split 0 

if vals [0] .upper 0 == "FRONT.NODES : " : 
n = int(vals[l]) 
while len (nodes) < n: 

buff = fdb.readlineO 

vals = buff. split 0 

for i in range (len (vals) ) : 

nodes . append (int (vals [i] ) ) 
f ront_nodes [f ] = nodes 
nodes = [] 
f+=l 

fdb . closeO 

if len(front_nodes) != nfronts: 

print "Only found Xi crack fronts in Xs" % (i+1 ,fdbFile) 

#print "\nHere is the front_nodes dict:\n" 

#print front_nodes 

# Second, create a nested dictionary where, for each crack front id (key), 

# the value is a dictionary with front node id’s as keys and corresponding 

# VecSD points located d distance "behind" each crack front node as values 

# 

# pts_per_front = {crack_front_id: {f ront_node_id: CTD_point}} 

# 

pts_per_front = {} 

# Loop over the crack fronts (handles multiple fronts) 
for front in front_nodes: 
pts = {} 

# For this crack front, loop over the crack front nodes 
for nid in front_nodes [front] : 
flag = 0 

adj_nodes = mesh. GetAdjacentCornerNodes (nid) 

# Rely on topology to find the adjacent node directly 

# "behind" the crack front node (this approach may break 

# down in certain instances, e.g. if crack front mesh 
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# template is not used.) The immediately adjacent 

# nodes on the main and mate surfaces should be initially 

# coincident, so just query the main surface node set for 

# a match. Note that for quadratic elements, only corner 

# nodes will return an adjacent corner node, 
for node in adj_nodes: 

if node in main_side_nodes : 

flag =1 # the front node is a corner node 

break 

if flag == 0: # Mid-side nodes along front don’t have 

continue # adjacent nodes on the crack face 

adj_ncoord = mesh. GetNodeCoords (node) 
front_ncoord = mesh. GetNodeCoords (nid) 

# Generate vector pointing from the crack front node 

# to the adjacent surface node 
vec = [] 

for i in range (3): 

vec . append (adj _ncoord [i] -f ront_ncoord [i] ) 
mag = sqrt(pow(vec[0] ,2)+pow(vec[l] ,2)+pow(vec[2] ,2)) 

# Now get point coordinates using the unit vector 

# components multiplied by the user specified distance, d. 

X = front_ncoord [0] + vec[0]/mag*d 

y = front_ncoord[l] + vec[l]/mag*d 
z = front_ncoord[2] + vec[2]/mag*d 
ptsEnid] = Vec3D. Vec3D(x,y ,z) 
pts_per_front [front] = pts 

print "\nReturning CTD(A) points at d=7og behind crack front nodes..." % d 
print pts_per_front 
print len(front_nodes) 
return pts_per_front 

def setTol(mesh,main_side_nodes) : 

# Set Pointinside tolerance to an adequate value local to the crack face: 

# 1/1000 of the smallest crack face edge length. This will be used in 

# the functions IsPointOutsideMesh (to check CTD pts) and GetElemsAtPoint 
min_length = mesh. GetMaxDimens ion () # A starting dimension 

print "\nSetting tolerance for MeshTools functions..." 
for nid in main_side_nodes : 

elengths = mesh. GetAdjacentSurfEdgeLengths (nid) 
min_length = min (min_length, min (elengths) ) 

tol = min_length/100 . 0 

mesh. SetPointlnsideTolerance(tol) 

return tol 

def getRelDispsAt InterestingNodes (disps_per_inc , nodes , j obName , curr_step , curr_inc , odb) : 

# Get displacements at the nodes of interest for the current load step/inc. 

# These are "relative" displacements because they do not account for initial 

# mapped displacements. 

print "\nAccessing odb to get displacements for nodes of interest" 
print "through load step 7oi increment 7oi..." 7o (curr_step , curr_inc) 

# Fill dictionary of [x,y,z] displacements for all relevant nodes 
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relative_disps = {} 


try: 

odbStep = odb . steps . keys 0 [curr_step-l] 
except : 

print "Looks like current step has not been written to odb file. Exiting." 
exit (0) 

# Loop through increments, from most recent stored increment to current 
for odbFrame in range (max (disps_per_inc) +1 , curr_inc+l) : 
try: 

fieldObject = odb . steps [odbStep] . frames [odbFrame] .f ieldOutputs [’U’] 
print "Increment 7oi..." 7o odbFrame 
for nid in nodes: 

values = fieldObject .values [nid- 1] . data 
relative_disps [nid] = values 
disps_per_inc [odbFrame] = relative_disps 
relative_disps = {} 
except : 

print "No odb info yet for Xs increment 7oi" 7o (odbStep, odbFrame) 

return disps_per_inc , odbStep 
print "Got ’em!" 


def getAbsDisps (mapped_disp ,relative_disp) : 

# Function to get absolute nodal displacements by summing the 

# mapped displacements with the current relative displacements. 

print "\nGetting absolute disps for nodes of interest at current step/inc... 
abs_disp = {} 

for nid in relative_disp : 
disp = [] 

for i in range (3): 

disp . append (relative_disp [nid] [i] +mapped_disp [nid] [i] ) 
abs_disp [nid] = disp 
return abs_disp 
print "Got ’em!" 

def updateDisplacements (mesh, abs_disp) : 

# Write a temporary file of absolute displacements for nodes of interest 

# to be read into mesh results (update to current displacements) 
print "\nUpdating current displacements..." 

tmpfile = open( ’tmp_current_disps . txt ’ , "w") 
for nid in abs_disp: 

line = ’7od\t7o. 15f\t7o. 15f\t7o. 15f \n’ 7o(nid,abs_disp[nid] [0] ,\ 

abs_disp [nid] [1] , abs_disp [nid] [2]) 

tmpfile. write (line) 
tmpfile . close 0 

mesh. ReadFeawdCpDisps ( ’ tmp_current_disps .txt ’ ) 
os . system( ’ del tmp_current_disps . txt ’ ) 

def computeAndWr iteCTDResult s (mesh , tol , d , CTDc , inpFile , f dbFile , j obName , odbStep , \ 
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curr_inc ,pts_per_front ,main_side_elems ,mate_side_elems) : 

# Compute CTD results and write to input file directory 
changeDirectory (inpFile) 

# Partition CTD values into subcritical, critical, and supercritical 

f lags_per_front = {} # A dictionary where key is crack front id 

# and value is list of flags for all CTD (A) points 
flag =0 # Flag written to results file 

# Write the file header 
resultsFile = ’CTDResults.txt’ 
if os .path. isf ile(resultsFile) : 

file = open(resultsFile , "a") 
else : 

file = open(resultsFile , "w") 
f ile.writeC ’CTD Results\n’) 

f ile. write ( ’Results provided along crack fronts (z_min to z_max)\n\n’) 

f Q 1 "t G ( ^ ^ 3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3|C3^ ^ ^ 

# Write sub-header 

file, write ( ’\n\nMesh file: 7os\n’ % inpFile) 
f ile. write ( ’Crack file: 7os\n’ % fdbFile) 
f ile. write ( ’Job name: 7os\n’ % jobName) 
file, write ( ’ABAQUS Step Name: 7os\n’ % odbStep) 
file, write (’Frame: 7oi\n’ % curr_inc) 

f ile. write ( ’Distance behind crack front nodes: 7og\n’ % d) 
f ile. write ( ’Critical opening displacement: 7og\n\n’ % CTDc) 

# Write CTD results for CTD points at each crack front 
for front in pts_per_f ront : 

print "front" 

flags = [] # A list of flags for each CTD (A) point 

file, write ( ’\nCrack Front ID: 7oi\n’ % front) 

file .write ( ’\nCorresponding\tCTD_I\t\tCTD_II\t\tCTD_III\tCTD_mag\tFlag’ ) 
file . write ( ’ \tCTD_ptZCoord\nFrontNodeId\n ’ ) 

f ile. write ( ’ \n’) 

temp_dict = {} 

print pts_per_f ront [front] 

# First sort the front node id’s by their z-coords. 

for nid in pts_per_f ront [front] : # Unsorted 

ncoords = mesh. GetNodeCoords (nid) 
zcoord = ncoords [2] 
temp_dict [zcoord] = nid 

sorted_front_nids = sortDictValsByKey (temp_dict) # Sorted 

# Move along crack front, computing CTD at points "behind" 

# each front node ID 

for nid in sorted_front_nids : 

CTD_pt = pts_per_front [front] [nid] 

status, dist = mesh. IsPointOutsideMesh(CTD_pt) 

if status != -2: 

file.write(’***CTD_pt C/og, 7og, 7og) is outside mesh and \ 

has been skipped!’ 7o (CTD_pt [0] ,CTD_pt [1] ,CTD_pt [2] ) ) 
f ile. write (’ (Behind front node 7oi)***’ 7o nid) 
continue 
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# Associate CTD_pt with elems on main/mate faces 

# (If less than two elems found, decrease meshTools tol and try again) 

# Reset MeshTools tol each time in case it was increased previously 
mesh . SetPointInsideTolerance (tol) 

elist = mesh.GetElemsAtPoint (CTD_pt) 
initial_tol = tol 

# Get displacement on the main side crack face for this CTD point 

# (displacement accounts for initial mapped displacements) 
main_eid = -1 

while main_eid == -1: 
for eid in elist: 

# Loop through elems until find main_side element 
if eid in main_side_elems : 

#print "main side elem for nid Xg: 7og" 7o (nid,eid) 

main_eid = eid 

break 

# If found main_side element, jump out of while loop 
if main_eid != -1: 

continue 

# Otherwise reset the tolerance and try again 
else : 

#print "newtol" 

new_tol = initial_tol*10 . 0 

#print new_tol 

mesh. SetPointInsideTolerance (new_tol) 
elist = mesh. GetElemsAtPoint (CTD_pt) 
initial_tol = new_tol 
displ = mesh.GetPtDisp(CTD_pt ,main_eid) 
ct=l 

# Get displacement on the mate side crack face for this CTD point 

# (displacement accounts for initial mapped displacements) 
mate_eid = -1 

while mate_eid == -1: 
for eid in elist: 

#print "elem:" 

#print eid 

# Loop through elems until find mate_side element 
if eid in mate_side_elems : 

#print "mate side elem for nid 7og: 7og" 7o (nid, eid) 

mate_eid = eid 

break 

# If found mate_side element, jump out of while loop 
if mate_eid != -1: 

continue 

# Otherwise reset the tolerance and try again 
else : 

#print "newtol" 

new_tol = initial_tol*10 . 0 

#print new_tol 

mesh. SetPointInsideTolerance (new_tol) 
elist = mesh. GetElemsAtPoint (CTD_pt) 
initial_tol = new_tol 
disp2 = mesh.GetPtDisp(CTD_pt ,mate_eid) 
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# Sometimes same element associated for main and mate surfaces 

# (notice especially with poorly shaped elements around crack front) . 

# Check for this condition and exit if that’s the case, 
if main_eid == mate_eid: 

print "Same main/mate elem found for CTD pt behind node Xg" 7o nid 
print "Check the quality of the mesh. Exiting!" 
exit (0) 

# Compute CTD’s and raise flag if critical value is reached 

# NOTE: The specific fracture criterion evaluation is handled in the 

# function evaluateFractureCriterion 
CTD = displ - disp2 # vector of disps 
CTD_I = abs(CTD[l]) # Mode I --> y-disp 
CTD.II = abs(CTD[0]) # Mode II--> x-disp 
CTD.III = abs(CTD[2]) # Mode III-> z-disp 

CTD.mag = sqrt (pow(CTD_I ,2)+pow(CTD_II ,2)+pow(CTD_III ,2) ) 
tol = CTDc*0.02 # Tolerance is 2% of specified critical value 

if abs (CTD_mag-CTDc) < tol: # CTD = CTDc +- tol 

flag = 1 

flags . append (flag) 

elif CTD.mag > CTDc: # CTD » CTDc 

flag = 2 

flags . append (flag) 
else : 

flag =0 # CTD « CTDc 

flags . append (flag) 

file.write(’yoi\t\tyog\tyog\tyog\tyog\t7oi\tyog\n’ 7o \ 

(nid , CTD_I , CTD_I I , CTD_I I I , CTD.mag , f lag , CTD.pt [2] ) ) 

flags_per_front [front] = flags 
print "\nDone writing to the results file!" 
file . close 0 
return f lags_per_f ront 

def evaluateFractureCriterion(f lags_per_f ront , eval) : 

# eval --> =0.0 corresponds to mid-thickness only option 

# --> !=0 corresponds to 7 of CTD points that must meet criterion 

frac_status_per_front = {} # A dictionary where key is crack front id 

# and value is: 

# 0 --> subcritical, continue 

# 1 --> critical, extend 

# 2 --> supercritical, back-up 

# For the mid-thickness only criterion evaluation: 
if eval == 0.0: 

# Loop through each crack front 
for front in f lags_per_front : 
index = [] 

# If even number of points through-thickness, 

# get indices of the middle two points 
if len(f lags_per_f ront [front] ) 7 2 == 0: 

index . append ( (len (f lags_per_f ront [front] ) -2) /2) 
index. append (len(flags_per_f ront [front] )/2) 
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# Otherwise, if odd number of points through-thickness, 

# get index of the middle CTD(A) point 
else : 

index . append ( (len(f lags_per_front [front] ) -1) /2) 
for i in index: 

# The critical case: 

if flags_per_front [front] [i] == 1: 
f rac_status_per_f ront [front] = 1 
break 

# The supercritical case: 

elif flags_per_f ront [front] [i] == 2: 
f rac_status_per_f ront [front] = 2 
break 

# The subcritical case: 

elif flags_per_f ront [front] [i] == 0: 
f rac_status_per_f ront [front] = 0 
break 

# For the evaluation case where evalX of points must meet criterion 
elif eval != 0.0: 

# Loop through each crack front 
for front in f lags_per_front : 

# Count number of sub-, super-, and critical flags along front 

sub =0 # Clear the count 

super =0 # Clear the count 

crit =0 # Clear the count 

total = len(f lags_per_f ront [front] ) 
for flag in f lags_per_f ront [front] : 
if flag == 1: 
crit+=l 

elif flag == 2: 
super+=l 
elif flag == 0: 
sub+=l 

# The critical case: 

if (float (crit) /float (total) ) >= eval: 
frac_status_per_f ront [front] = 1 
continue 

# The supercritical case: 

elif (float (super) /float (total) ) >= 1-eval: 
frac_status_per_f ront [front] = 2 
continue 

# The subcritical case: 
else : 

frac_status_per_f ront [front] = 0 


return f rac_status_per_front 

# NOTE: This script script cannot be executed from the CAE if *MAP SOLUTION is 

# specified since the keyword is not yet supported by the CAE (and consequently 

# the old job cannot be specified) 

if name == " main ": 
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print doc 


inpFile = sys.argv[-4] # deformed . inp file (current crack step) 

undefInpFile = sys.argv[-3] # undeformed .inp file in Undeformed sub-dir 

fdbFile = sys.argv[-2] # .fdb file generated by FRANC3D\NG 

odbFile = sys.argv[-l] # .odb file of previous analysis that ABAQUS 

# should use to map solution 

changeDirectory (inpFile) 

if not odbFile .upper 0 == "NONE": 

# Old job path 

oldJob = odbFile.splitC'.") [0] 
mapSoln = "Y" 

# Check for all the necessary restart files of the old job 
res = oldJob+ ’ . res ’ 

mdl = oldJobH-’ .mdl’ 
stt = oldJobH-’ . stt ’ 
prt = oldJob+’ .prt ’ 
if not os .path. exists (res) : 

print "\tCannot find restart file for the oldJob." 
print "\tCheck that it exists. Exiting." 
exit (0) 

if not os .path. exists (mdl) : 

print "\tCannot find .mdl file for the oldJob." 
print "\tCheck that it exists. Exiting." 
exit (0) 

if not os .path. exists (stt) : 

print "\tCannot find stt file for the oldJob." 
print "\tCheck that it exists. Exiting." 
exit (0) 

if not os .path. exists (prt) : 

print "\tCannot find restart file for the oldJob." 
print "\tCheck that it exists. Exiting." 
exit (0) 

else : 

oldJob = "None" 
mapSoln = "N" 

# Check that status file doesn’t already exist. If it does, delete it so that 

# we don’t tail the existing file 
staFile = inpFile . split (".") [0] +". sta" 
if os .path. isf ile(staFile) : 

response = raw_ input ("\nThe file Xs already exists. Overwrite file? (y/n) : "\ 
7o StaFile) 

if response .upper 0 [0] == "N" : 
print "Exiting program." 
exit (0) 

elif response .upper 0 [0] == "Y" : 
os . system( ’del ’+staFile) 
else : 

print "Unknown response. Exiting." 
exit (0) 
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# Request CTD parameters: 

print " USER INPUTS " 

# --> Distance, d, behind crack front 

d = askForFloat ("\nMonitor opening at this distance behind crack front: ") 

# --> CTD or CTOA 

response = raw_ input ("\nWould you like to specify critical ’CTD’ or ’CTOA’?: ") 
while response .upper 0 not in ("CTD" , "CTOA") : 

response = raw_input ("\n\tSpecify ’CTD’ or ’CTOA’: ") 
if response .upper 0 in "CTD": 

CTDc = askForFloat ("\n\tCritical CTD: ") 
elif response .upper 0 in "CTOA": 

CTOAc = askForFloat ("\n\tCritical CTOA (degrees): ") 

# Convert CTOAc value to CTDc for Mode I 
rad = CT0Ac*math.pi/180 
CTDc = rad*d 

# --> Criterion 

print "\nHow is fracture criterion evaluated?" 
response = askYesNo("\n\tEvaluate at mid-thickness only?: ") 
if response in "Y" : 
eval = 0.0 
if response in "N" : 

print "\n\tCriterion will be evaluated at a percentage of points" 

response = askForInt ("\n\tSpecify percentage of points that must meet criterion: ") 
eval = float (response) /lOO 

print "\n STARTING FEA " 

# Start the analysis 

runAnalysis (d , CTDc , eval , oldJob ,mapSoln , undef InpF ile , inpF ile , f dbFile) 


B. Deformation Mapping Script 

Two scripts are provided for mapping deformation from the previous mesh to the current, 
undeformed mesh generated by FRANC3D\NG. The first script, called runMapScript.py^ simply 
prompts the user for a series of input, including paths to various required files. The user will also 
be asked whether or not the analysis is a local-global analysis (i.e. if the model has global and 
sub-regions). If it is, then displacements on the global region will simply be transferred from old 
to new meshes, as the global mesh remains unchanged and does not require invoking the inverse 
isoparametric mapping routine. The user will also be asked whether or not the previous analysis is 
a restart analysis so that the correct restart files are accessed and displacements are taken from the 
correct load increment. The script can be run from any directory. The user should execute the script 
from the MS-DOS command window using a command like: “...>abaqus python runMapScript.py”. 

The script runMapScript.py automatically calls the script master Map Displacements, py^ which 
retrieves displacements from the previous FE analysis, maps them onto the current undeformed 
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mesh, and generates a deformed mesh file (*.inp) for the current crack increment. The deformed 


*.inp file will contain the *Map Solution ABAQUS [1] keyword and an equilibration step, during 
which all nodes with applied boundary conditions are fixed. The script also generates the sub- 
directory Undeformed and moves the original undeformed mesh into that directory. Output and 
print statements during mapping are written to the abaqus.rpy file in the current working directory. 

After mapping, the script CTDjohControlCalLpy is called and the next FE analysis automati- 
cally begins and continues until CTD^CTDcrit- 

Successful execution of these scripts require the compiled libraries Vec3D, ColTensor, Full- 
Tensor, and MeshTools to be located in a directory accessible by ABAQUS® (e.g. in the 
ABAQUS >Python>Obj directory). Elements currently supported for mapping include linear and 
quadratic brick elements (C3D8 and C3D20), quadratic wedge elements (C3D15), and quadratic 
tetradhedral elements (C3D10). 

Electronic versions of the deformation mapping scripts are available for download at 
WWW. efg. eornell. edu. 


? ? ? 

runMapScript .py 
This script: 

1) generates a .model file from the current undeformed 
. inp file, 

2) generates mapped displacements from old mesh to new 
undeformed mesh, 

3) appends these displacements to the .model file, 

4) moves the undeformed .inp file to a new directory, and 

5) rewrites the .inp file in the deformed configuration, 

6) automatically calls CTDjobControlCall .py to start the 
nextFEA if "yes" at prompt. 

Requires an initial undeformed .inp file for the current 
step, a .model file for the previous step (with a header 
for mapped displacements) , as well as the following compiled 
scripts : 

Vec3D, ColTensor, FullTensor, MeshTools 

AD Spear, MG Veilleux, and JD Hochhalter 

(written October 2009 with subsequent modifications) 

? ? ? 

import os , sys 
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from sys import exit 

from abaqusConstants import * 

import odbAccess 

def askForPath (prompt , extension) : 

# Request path from user 
response = raw_ input (prompt) 
while not os .path. exists (response) : 

print "Specified path could not be located." 
print "Check that path exists and try again." 
response = raw_ input (prompt) 
while not response . split (’.’) [-1] == extension: 

print "\n**File path should have extension ’7oS’**" % extension 
response = raw_ input (prompt) 
return response 

def askYesNo (question) : 

# Request yes/no response from user 

# (Returns "Y" or "N") 
response = raw_ input (question) 

while response .upper 0 [0] not in ("Y","N"): 
response = raw_ input (question) 
if response .upper 0 [0] not in ("Y","N"): 
print "Wrong answer! Try again" 
return response .upper () [0] 

def getFdb(path) : 
fdb = str() 

# First try splitting on "." 
tmp = path. split (". ") 
tmp[-l] = ’.fdb’ 

for i in range (len (tmp) ) : 
f db+=tmp [i] 

# If doesn’t exist, try splitting on "_" 
if not os .path. exists(fdb) : 

fdb = str() 
tmp = path. split ("_") 
tmp[-l] = ’.fdb’ 
for i in range (len (tmp) ) : 
f db+=tmp [i] 

if not os .path. exists (fdb) : 

print "The file Xs cannot be located." % fdb 
exit (0) 

print "\nGetting Xs" % fdb 
return fdb 

def getOdb(path, localglobal ,res) : 
odb = str() 
tmp = path. split (". ") 
if localglobal == "Y" and res == "N" : 
tmp[-l] = ’_full.odb’ 
elif localglobal == "Y" and res == "Y" : 
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tmp[-l] = ’ _full_res . odb’ 
elif localglobal == "N" and res == "Y" : 
tmp[-l] = ’_res.odb’ 
else : 

tmp[-l] = ’.odb’ 
for i in range ( len (tmp) ) : 
odb+=tmp [i] 

while not os .path. exists(odb) : 

print "The file Xs cannot be located." % odb 
odb = raw_ input ("\nEnter path to corresponding odb file:\n") 
print "\nGetting Xs" % odb 
return odb 

def askForInt (number) : 

# Request int from user 
response = int (raw_ input (number) ) 
while response < 0 : 

print "Requested value must be a positive integer. Try again." 
response = int (raw_input (number) ) 
return response 

if name == " main ": 

print doc 

currInpFile = "" 
currInpFile_full = "" 

# Ask user if analysis is local-global analysis 
promptO = "Is analysis local-global \ 

(i.e. *_full files generated by F3DNG)?:\n" 
localglobal = askYesNo (promptO) 

# Request . inp file of job from which displacements will be mapped 

# (if local-global, just the local input file generated by F3DNG) 
if localglobal == "Y" : 

prompt 1 = "\nPath of old, undeformed input file \ 

(local region) \n(e .g. C : / . . . /oldMesh. inp) : \n" 
previnp = askForPath(promptl , ’ inp ’ ) 
else : 

prompt 1 = "\nPath of old, undeformed input file \ 

\n(e.g. C: /... /oldMesh. inp) : \n" 
previnp = askForPath(promptl , ’ inp ’ ) 

# Get top level directory to search for .fdb and .odb files 
tmp = previnp . split ("/Undeformed/") 

dir = tmp [0] +"/"+tmp [-1] 
print dir 

# Search for corresponding .fdb file 
fdbPrev = getFdb(dir) 

print fdbPrev 
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# If the analysis is local-global, need both 

# the * . inp and *_full.inp files 
if localglobal == "N" : 

# Request current .inp file to which displacements will be mapped 
prompt2 = "\nPath of current .inp file to which displacements \ 
will be mapped\n(e . g. C: / . . . /newMesh. inp) : \n" 
currInpFile = askForPath(prompt2, ’ inp’ ) 

# Check for existence of corresponding .fdb file 
fdbCurr = get Fdb (currInpFile) 

elif localglobal == "Y" : 

# Request * . inp generated by F3DNG 

prompt2a = "\nPath of most recent (local) * . inp file \ 
generated from F3DNG:\n" 
currInpFile = askForPath(prompt2a, ’ inp ’ ) 

# Check for existence of corresponding .fdb file 
fdbCurr = get Fdb (currInpFile) 

# Also check for the *_full.inp of the current, joined file 
tmp = currInpFile . split (". ") 

tmp[-l] = "_full.inp" 
for i in range (len (tmp) ) : 

currInpFile_full+=tmp [i] 
if not os .path. exists (currInpFile_full) : 

print "Cannot locate the file 7oS" % currInpFile_full 
exit (0) 
else : 

print "Retrieving Xs" % currInpFile_full 

# Search for .odb file of old mesh and access the file 
promptO = "\nMap disps. from a *_full_res . odb file?:\n" 
res = askYesNo (promptO) 

odbFile = getOdb (dir , localglobal , res) 
try: 

odb = odbAccess . openOdb(path=odbFile ,readOnly=TRUE) 
except : 

print "ERROR: Unable to open the specified odb Xs" % odbFile 
exit (0) 

# Ask for analysis step to map displacements from 

prompt3 = "\nName of step displacements will be mapped from:\n" 
stepName = raw_ input (prompt 3) 
if StepName not in odb . steps .keys () : 

print "\nERR0R: Step ’Xs’ does not exist in 7os\n"\ 

"\tCheck for the case in the step name." % (stepName, odbFile) 
exit (0) 

# Ask for frame number (load increment) to map displacements from 
frame = askForInt ("\nFrame number (increment) to map \ 

displacements from:\n") 
if len (odb . steps [stepName] . frames) < frame: 

print "\nERR0R: Frame Xi does not exist in 7os\n"\ 

"\tCheck for the case in the .sta file." % (frame, odbFile) 
exit (0) 

odb . closeO 


64 



# Input parameters to masterMapDisplacements .py script*** 
os . system ( ’abaqus cae noGUI=masterMapDisplacements .py -- ’\ 

+res+’ ’+localglobal+’ ’+prevlnp+’ ’+fdbPrev+’ ’+odbFile+’ ’\ 
+currInpFile+’ ’+currInpFile_full+ ’ ’+fdbCurr+’ ’\ 

+stepName+’ ’+str (frame) ) 

# Pass arguments to CTDjobControlCall script 
prompt4 = "\nStart FEA with CTD job control? :\n" 
start = askYesNo(prompt4) 

if start == "Y" : 

# Current deformed input file 
if localglobal == "Y" : 

defInpFile = currInpFile_full 
else : 

defInpFile = currInpFile 

# Undeformed input file located in Undeformed sub-dir 
undefInpFile = str() 

tmp = currInpFile . split ("/") 
fileName = tmp[-l] 
tmp[-l] = "Undeformed" 
for i in range (len (tmp) ) : 

undef InpFile+=tmp [i] +"/" 
undef InpFile+=f ileName 

# .fdb file 

fdb = getFdb( currInpFile) 

print "Calling CTDjobControlCall .py" 

os . system (’ abaqus python CTDjobControlCall .py -- ’\ 

+def InpFile+ ’ ’+undef InpFile+’ ’+fdb+’ ’+odbFile) 

else : 

print "Displacement mapping complete!" 

# Finally, move the abaqus.rpy file to the directory of the current mesh 
print "Moving the abaqus.rpy file..." 

currPath = str() 
tmp = currInpFile . split ("/") 
for i in range (len (tmp) -1) : 
currPath+=tmp [i] 
currPath+="/" 
currPath+="abaqus . rpy" 
os . rename ( "abaqus . rpy" , currPath) 
os .remove ("abaqus_ac is . log") 


? ? ? 

masterMapDisplacements .py 
Executed by runMapDisplacements .py 
? ? ? 

import Vec3D 
import ColTensor 
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import FullTensor 

import MeshTools 

import os, glob, sys, time 

from abaqus import * 

from abaqusConstants import * 

import part 

import assembly 

import mesh 

from sys import argv, exit 
import odbAccess 

def changeDirectory (pathToFile) : 
path = str() 

tmp = pathToFile . split ("/") 
for i in range (len (tmp) -1) : 
path+=tmp [i] 
path+="/" 
print path 
os . chdir (path) 

def createSubDirectory (f ilePath, subDirName) : 
newdir = ’ ’ 

path = filePath. split ("/") 
for i in range (len (path) -1) : 
newdir+=path[i] +"/" 
newdir+=subDirName+"/" 
if not os .path. isdir (newdir) : 
os . mkdir (newdir) 
return newdir 

def createMeshToolsObject (filePath, extension) : 
base = str() 

tmp = filePath. split (". ") 
tmp [-1] = ’ ’ 

for i in range (len (tmp) ) : 
base+=tmp [i] 

meshObj = MeshTools . MeshTools (base , extension) 
return meshObj 

def getGroupList (meshObj , listName , type) : 

list = meshObj .GetGroupInfo (listName, type) 
if len(list) == 0: 

print "***No members found in group: Xs ***" % listName 
return list 

def printTimeStamp(start_time) : 

# Print time stamp 

end_time = time.timeO 

mins = int ( (end_time-start_time) /60) 

secs = round( (end_time-start_time)7o60,3) 
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print "\nTotal time for mapping: 7oi:7og" 7o (mins, secs) 

# ******************* get initial displacements from old mesh ********************* 

def getInitialDisplacements (f ilePath) : 

# Get initial displacements, if any, from the ## Displacement header of the 

# .model file of the previous mesh 
initial_disps = {} 

disp_file = str() 
buff = filePath. split ( V’ ) 
for i in range(len(buf f ) -1) : 
disp_f ile+=buf f [i] + V ’ 
disp_f ile+= ’mapped_disp . txt ’ 
print disp_file 

# If the disp.txt file doesn’t exist, continue, but inform analyst 
if not os .path. exists(disp_f ile) : 

print "***Could not locate mapped_disp.txt from previous analysis!***" 
print "(Make sure the file is located with the previous mesh file)" 
return initial_disps 

file = open(disp_f ile, ’r ’ ) 
buff = file. readline 0 
while buff : 

data = buff. split 0 
if data [0] .upper 0 [0 : 4] == "DISP" : 
buff = file .readline 0 
data = buff. split 0 
while buff : 

nid = int(data[0]) 
tmp = [] 

for i in range (3): 

tmp . append (float (data[i+l] ) ) 
initial_disps [nid] = tmp 
buff = file. readline 0 
data = buff. split 0 
buff = file .readline 0 

file. closeO 

print "***Retrieved initial displacements for 7od nodes***" 7o len(initial_disps) 
return initial_disps 

# **************** get relative displacements from previous analysis ************** 

def getRelativeDisplacements (odbFile , stepName , frame , oldNodeList , localglobal) : 

# Extract displacements at the critical step and increment from the previous 

# analysis 

relative_disps = {} # Dictionary containing odb disps. for all nodes 

l_relative_disps = {} # Dictionary containing odb disps. for local nodes 

g_relative_disps = {} # Dictionary containing odb disps. for global nodes 

fullNodeList = [] # List of all nodes from odb-->indices correspond to 

# indices in .odb node list 

g_odbIndices = [] # List of node indices from odb global nodes for 

# referencing nodes by index rather than node ID 

# (if local-global analysis) 
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try: 

odb = odbAccess . openOdb(path=odbFile ,readOnly=TRUE) 
except : 

print "ERROR: Unable to open 7oS" % odbFile 

exit (0) 

# Get step number corresponding to stepName 
for i in range (len (odb . steps .keys ())) : 

if StepName in odb . steps . keys () [i] : 
stepNum = i+1 

# Create list of all nodes in model 

# NOTE: nodeLabels do not necessarily correspond to fieldObject index 

# especially for local/global join 

for instance in odb . root Assembly . instances . keys ( ) : 

for i in range(len(odb.rootAssembly. instances [instance] .nodes) ) : 

f ullNodeList . append (odb . root Assembly . instances [instance] . nodes [i] . label) 

fieldObject = odb. steps [stepName] .frames [frame] .fieldOutputs [’U’] 

# If the model is local-global, then separate odb displacements into 

# local and global dictionaries, to be handled separately 
if localglobal == "Y" : 

for nid in fullNodeList : 

index = fullNodeList . index (nid) 
disps = fieldObject .values [index] . data 
if nid in oldNodeList: 

l_relative_disps [nid] = disps 
else : 

g_relative_disps [nid] = disps 
g_odbIndices . append (index) 

# Otherwise include all odb disps into a single dictionary 
elif localglobal == "N" : 

for nid in fullNodeList: 

index = fullNodeList . index (nid) 
disps = fieldObject .values [index] . data 
relative_disps [nid] = disps 


else : 

print "Unknown response to local/global command prompt" 
exit (0) 

odb . closeO 

if len(relative_disps) > 0: 

print "***Retrieved odb disps for Xd nodes" % len(relative_disps) 
elif len(l_relative_disps) > 0 and len(g_relative_disps) >0 : 

print "***Retrieved odb disps for Xd local nodes and Xd global nodes" % \ 
(len(l_relative_disps) , len(g_relative_disps) ) 

else : 

print "***No relative displacements retrieved from odb. \ 

Exiting map script." 

exit (0) 

return stepNum, fullNodeList ,relative_disps , l_relative_disps , \ 
g_relative_disps ,g_odbIndices 
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# ************************ COMPUTE ABSOLUTE DISPLACEMENTS ************************* 

def computeAbsoluteDisplacements (initial_disps ,relative_disps ,nodeList ,prevMesh) : 

# Compute the absolute nodal displacements for the old mesh at the critical 

# load step/ increment of interest 
abs_disps = {} 

# Quick check that length of lists are the same 
if len(relative_disps) != len(nodeList) : 

print "Size of lists disagree in computeAbsoluteDisplacements function!" 
exit (0) 

if len(initial_disps) != len(relative_disps) : 
if len(initial_disps) == 0: 

print "***No initial displacements found in previous mesh***" 
print "***Check that previous mesh did not have initial \ 
displacements***" 

print "*** (Mapping will only include relative displacements from the odb) " 

# Set initial displacements to zero 
for node in nodeList: 

initial_disps [node] = [0,0,0] 

else : 

print "***Initial and relative displacements found for different \ 
number of nodes!" 

print "***Check model and odb files from previous mesh for the discrepancy!" 
exit (0) 

for node in nodeList: 
d = [] 

for i in range (0,3): 
try: 

val = float (initial_disps [node] [i] )+f loat (relative_disps [node] [i] ) 
d. append (val) 
except : 

print "***Error computing absolute displacements for node 7od!" % node 
print "***Check that initial and relative disp vals exist for the node" 
exit (0) 

abs_disps [node] = d 
return abs_disps 

# * * * * * * * * * * * * * * * * * * * map displacements to new undeformed mesh ******************** 

def mapDisplacementsToUndef ormedMesh(meshl ,abs_disps ,mesh2) : 

# sic************** mapDisplacementsToUndef ormedMesh SUB-FUNCTIONS ************** 

? ? ? 

def par seModelFileForGroup (file, keyword) : 
list = [] 

fd = open(f ile , "r") 

buff = fd. readline 0 
while buff : 
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data = buff. split 0 

if len(data) == 5: # Implying group sub-header line 

if data [4] . upper () == keyword: 

print "Found Xs in 7oS" % (keyword, file) 
buff = fd. readline 0 
data = buff. split 0 
while buff and data[0] != 

for i in range (len (data) ) : 

list . append (int (data[i] ) ) 
buff = fd. readline 0 
data = buff. split 0 
break 

buff = fd. readline 0 
f d. close 0 
return list 

? ? ? 

def makeElemSet (mesh,nodeList) : 

# Create a set of elements that contain the nodes in nodeList 
setName = set() 

for nid in nodeList: 

elist = mesh. Get AdjacentElems (nid) 
for eid in elist: 

setName . add(eid) 

if len(setName) == 0: 

print "***ERR0R: None found for the case Xs" % set 
return setName 

def setPointInsideTol(mesh,mainNodeList ,mateNodeList) : 

# Set point inside tolerance to 1/1000 of the smallest crack 

# face edge length 

min_length = mesh. GetMaxDimens ion () # Initialize a minimum length 

for nid in mainNodeList : 

elengths = mesh. Get AdjacentSurfEdgeLengths (nid) # function returns 
min_length = min (min_length, min (elengths) ) # error if nid is 

# not surface node 

for nid in mateNodeList : 

elengths = mesh. Get AdjacentSurfEdgeLengths (nid) 
min_length = min (min_length, min (elengths) ) 

print "min_length" 

print min_length 

initial_tol = min_length/1000 . 0 

mesh. SetPointInsideTolerance(initial_tol) 

print "***Min edge length is: /of" 7o min_length 

print "***Tolerance is: /of" 7o initial_tol 

return initial_tol 
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def mapCrackFaceNodeDisps 0 : 


print "***WARNING: main- and mate-side choices are assumed to be\n \ 
consistent in both crack models. " 
print "***Check fdb files to check validity of this assumption!" 

# The mesh objects for the previous and current meshes should contain 

# the necessary main/mate side node groups (excluding crack front nodes) 

# for the cases with and without a crack front template. 

# First collect the main/mate side node sets 
try: 

main_side_nodesl = getGroupList (meshl , "main_side_nodes" , "node") 
mate_side_nodesl = getGroupList (meshl , "mate_side_nodes" , "node") 
except : 
try: 

main_side_nodesl = getGroupList (meshl , "all_main_side_nodes" , "node") 
mate_side_nodesl = getGroupList (meshl , "all_mate_side_nodes" , "node") 
except : 

print "***No main/mate side nodes found for previous mesh!***" 
exit (0) 


try: 

main_side_nodes2 = getGroupList (mesh2 , "main_side_nodes" , "node") 
mate_side_nodes2 = getGroupList (mesh2 , "mate_side_nodes" , "node") 
except : 
try: 

main_side_nodes2 = getGroupList (mesh2 , "all_main_side_nodes" , "node") 
mate_side_nodes2 = getGroupList (mesh2 , "all_mate_side_nodes" , "node") 
except : 

print "***No main side nodes found for current mesh!***" 
exit (0) 

#print len(main_side_nodesl) 

#print len(mate_side_nodesl) 

#print len(main_side_nodes2) 

#print len(mate_side_nodes2) 

# Next make sets of main/mate crack face elements for the previous mesh 
print "***Collecting main/mate side elems from previous mesh..." 
main_side_elemsl = makeElemSet (meshl ,main_side_nodesl) 
mate_side_elemsl = makeElemSet (meshl ,mate_side_nodesl) 

# Set an initial Pointinside tolerance based on the previous mesh 

# This is the tolerance for checking existence of a point from 

# the new mesh in the old mesh. 

initial_tol = setPointInsideTol (meshl ,main_side_nodesl ,mate_side_nodesl) 

# Map displacements onto the crack face nodes of the current mesh 
crkf ace_node_disps = {} 

# First map displacements to main_side nodes of current mesh 
for nid in main_side_nodes2 : 

# Reset tolerance in case it was previously increased, 
tol = initial_tol 

meshl . SetPointInsideTolerance (tol) 
elist = [] 
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info = mesh2 .GetNodelnfo(nid) 
pt = info[0] 

status, dist = meshl . IsPointOutsideMesh(pt) 

# Make sure pt . is within the bounding box of the mesh 
if status == -3: 

print "***Node id: %i " 1 nid 
print "***Nodal coords: " 
print pt 

print "***Status (old eid) : 7oi" 7o status 

raise ValueError, "Query point must be inside mesh" 

# If the point in the new mesh exists in a crackface element of 

# the old mesh, then interpolate displacements only from the containing 

# element (i.e. do not interpolate over surrounding elements in high 

# gradient region!) 

# Associate node with elements in previous mesh. If empty range tree 

# search (which happens numerically on rare instances) , increase tol 

# and try again, 
while len(elist) == 0: 

try: 

elist = meshl .GetElemsAtPoint(pt) 
except : 

print "Tol increased for nid Xg due to empty range tree search" \ 
7o nid 

new_tol = tol*10.0 

meshl . SetPointInsideTolerance (new_tol) 
tol = new_tol 

main_eid = -1 
mate_flag = 0 

# Loop through associated elems and look for containing main_side elem 

# in previous mesh (won’t be the case for nodes in new surface area) . 
for eid in elist: 

if eid in main_side_elemsl : 
main_eid = eid 
break 

if eid in mate_side_elemsl : 
mate_flag = 1 

# If a mate_side element was associated, but not a main_side element, 

# then enter while loop to increase query tolerance, 
while main_eid == -1 and mate_flag == 1: 

new_tol = tol*10.0 

print "increasing tolerance to Xf for nid Xg" 7o (new_tol ,nid) 
meshl . SetPointInsideTolerance (new_tol) 
elist = meshl . GetElemsAtPoint (pt) 
for eid in elist: 

if eid in main_side_elemsl : 
main_eid = eid 
break 
tol = new_tol 
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disp = meshl .GetPtDispCpt ,main_eid) 
crkf ace_node_disps [nid] = disp 
print "looped through all main_side_nodes2" 

# Next map displacements to mate_side nodes of current mesh 

# (just like above) 

for nid in mate_side_nodes2 : 

# Reset tolerance in case it was previously increased, 
tol = initial_tol 

meshl . SetPointInsideTolerance (tol) 
elist = [] 

info = mesh2 .GetNodelnfo(nid) 
pt = info[0] 

status, dist = meshl . IsPointOutsideMesh(pt) 

# Make sure node clearly exists inside previous mesh 
if status == -3: 

print "***Node id: %i " 1 nid 
print "***Nodal coords: " 
print pt 

print "***Status (old eid) : 7oi" 7o status 

raise ValueError, "Query point must be inside mesh" 

# Associate node with elements in previous mesh. If empty range tree 

# search (which happens numerically on rare instances) , increase tol 

# and try again, 
while len(elist) == 0: 

try: 

elist = meshl .GetElemsAtPoint(pt) 
except : 

print "Tol increased for nid Xg due to empty range tree search" \ 
7o nid 

new_tol = tol*10.0 

meshl . SetPointInsideTolerance (new_tol) 
tol = new_tol 

mate_eid = -1 
main_flag = 0 

# Loop through associated elems and look for containing mate_side elem 

# in previous mesh (won’t be the case for nodes in new surface area) . 
for eid in elist: 

if eid in mate_side_elemsl : 
mate_eid = eid 
break 

if eid in main_side_elemsl : 
main_flag = 1 

# If a main_side element was associated, but not a mate_side element, 

# then enter while loop to increase query tolerance, 
while mate_eid == -1 and main_flag == 1: 

new_tol = tol*10.0 

print "increasing tolerance to Xf for nid Xg" 7o (new_tol ,nid) 
meshl . SetPointInsideTolerance (new_tol) 
elist = meshl . GetElemsAtPoint (pt) 
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for eid in elist: 

if eid in mate_side_elemsl : 
mate_eid = eid 
break 
tol = new_tol 

disp = meshl .GetPtDispCpt ,mate_eid) # No averaging over adjacent elems 
crkf ace_node_disps [nid] = disp 
print "looped through all mate_side_nodes2" 
return crkf ace_node_disps , initial_tol 


# * * * * * * * * * * * * * * * * mapDisplacementsToUndef ormedMesh MAIN ******************* 
mapped_disps = {} 

# Note: Displacements currently associated with meshl are initial disps, 

# which were mapped during previous run of this script, 
for nid in abs_disps: 

coord = Vec3D . Vec3D(abs_disps [nid] [0] , abs_disps [nid] [1] , abs_disps [nid] [2]) 
meshl . UpdateNodalDisps (nid , coord) 

# Map node displacements at the crack face nodes first 

# (this routine interpolates displacements for points only within the 

# containing crackface element and doesn’t average over surrounding elements) 
(crkf ace_node_disps , initial_tol) = mapCrackFaceNodeDisps () 

# Now map node displacements for all nodes 
mesh2nodes = mesh2 . GetNodeList () 
mesh2nodes . sort () 

# Note that pt_query tolerance may be increased if a node is not deemed 

# inside the previous mesh; if tolerance is increased to the min. edge 

# length, an error is raised and the script is terminated, 
for nid in mesh2nodes: 

tol = initial_tol 

if nid in crkf ace_node_disps : 

mapped_disps [nid] = crkf ace_node_disps [nid] 
else : 

info = mesh2 .GetNodelnfo(nid) 
pt = info[0] 

status, dist = meshl . IsPointOutsideMesh(pt) 

# Status: -2 --> query _pt is clearly inside previous mesh 

# -1 --> query_pt is close to surface, but not clearly inside 

# -3 --> query_pt is not inside previous mesh 

while status != -2: 

if status == -1: 

new_tol = tol*10.0 

# (Uncomment these lines to put a limit on search tol) 

#if new_tol > initial_tol*1000 . 0 : 

#print "MeshTools tolerance exceeded smallest edge length!" 
#print "Node id Xi of currMesh is within Xf of prevMesh \ 
surface" % (nid, dist) 

#exit (0) 

meshl . SetPointInsideTolerance (new_tol) 
status, dist = meshl . IsPointOutsideMesh(pt) 
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tol = new_tol 

print "Tolerance has been reset to Xf due to nid Xi of currMesh"\ 
7o (tol, nid) 

else : 

txt = "Query point must be inside mesh" 
raise ValueError, txt 

# If the tolerance was increased for nid, reset to original 
meshl . SetPointInsideTolerance (initial_tol) 

# Get displacement by interpolating displacements from prevMesh 

# at point of interest 

mapped_disps [nid] = meshl . GetPtDisp(pt) 
return mapped_disps 


# * * * * * * * * * * * * * * * * * * APPEND MAPPED DISP RESULTS TO .MODEL FILE ******************** 

def writeMappedDispsToFile (mapped_disps , currInpFile , stepName , frame) : 

# Write the mapped displacement info to a .txt file 

# Create mapped_disp.txt in same directory as current input file 
changeDirectory (currInpFile) 

file = open ("mapped_disp. txt" , "w") 

file. write ("Displacements mapped from: Xs Increment: 7oi\n" % (stepName , frame) ) 
for nid in mapped_disps : 

f ile.write("7oi 7o.l6e 7o.l6e 7o.l6e\n" % (nid,mapped_disps [nid] [0] , \ 

mapped_disps [nid] [1] ,\ 
mapped_disps [nid] [2] ) ) 

file. closeO 

# * * * * * * * * * * * * * * * * * * * * * * * COMPUTE DEFORMED COORDINATES ************************** 

# * * * * * * * * * * * * * * * * For Global Region (if applicable) ************** 

def computeGlobalCoords (odbFile , stepName , frame , g_relat ive_disps , g_odbIndices) : 

# Create dictionary of final coordinates for global nodes in previous mesh 

# g_new_coords = {nid: (x,y,z)} 

print "Computing absolute disps for global region" 
g_new_coords = {} 

if len(g_relative_disps) != len(g_odbIndices) : 

print "Length of global node index list and rel_disp diet inconsistent!" 
eixt (0) 


try: 

odb = odbAccess . openOdb(path=odbFile ,readOnly=TRUE) 
except : 

print "ERROR: Unable to open Xs" % odbFile 
exit (0) 

# Initiate a counter for g_new_coords diet 
ent = 1 

# Simply add relative disps. to node coords to get deformed coords 

# for global region of previous mesh (no inverse isoparametric 

# mapping involved) 
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for instance in odb . root Assembly . instances . keys ( ) : 
instance = odb.rootAssembly. instances [instance] 
for index in g_odbIndices : 
tmp = [] 

node = instance .nodes [index] 
nid = node. label 
for i in range (0,3): 

val = float (g_relative_disps [nid] [i] )+float (node . coordinates [i] ) 
tmp . append (val) 

g_new_coords [cnt] = tmp 
cnt+=l 
odb . closeO 

return g_new_coords 

# * * * * * * * * * * * * * * * * For Local or Entire Region (whichever applicable) ************** 
def computeDeformedCoords(currInpFile,mapped_disps) : 

# Change directories to location of current . inp file if not there already 
changeDirectory (currInpFile) 

undef_coords = {} 
def_coords = {} 
scale = 1 

# Parse the original .inp for the current mesh to get the undeformed coords, 
file = open(currInpFile , "r") 

line = file. readline 0 
while line: 

if line [0 : 5] .upper 0 == "*N0DE": 
line = file .readline 0 

while line [0:1] != 

data = line . split (", ") 
nid = int(data[0]) 

undef _coords [nid] = [float (data[l] ) , float (data [2] ) , float (data [3] )] 
line = file . readline 0 
break 

line = file .readline 0 
file . close 0 

# Do a quick check that min/max node id’s are consistent 
if max (undef _coords) != max(mapped_disps) or \ 

min (undef _coords) != min(mapped_disps) : 

print "***Node numbering is inconsistent between .inp file and \ 
mapped_disps dict!\n" 

exit (0) 

# Compute deformed nodal coordinates 
for nid in undef _coords : 

x_new = undef _coords [nid] [0] +mapped_disps [nid] [0] *scale 
y_new = undef _coords [nid] [1] +mapped_disps [nid] [l]*scale 
z_new = undef _coords [nid] [2] +mapped_disps [nid] [2] *scale 
def _coords [nid] = [x_new,y_new,z_new] 
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return def_coords 


# * * * * * * * * * * * * * * * * * * * * * * * * * * WRITE DEFORMED MESH ( . inp) *************************** 
def wr iteDef ormedInpFile (def _coords , g_def _coords , currInpF ile , stepNum, frame , res) : 

# *************** Scan input file for BC nodesets SUB-FUNCTION *************** 
def getBCnodesets (inpFile) : 

# Scan for BC nodesets 
BCnodes = [] 

if ile = open(inpFile , "r") 
line = if ile. readline 0 

while line [0 : 6] .upper 0 != "*B0UND": 
line = if ile .readline 0 

while line: 

line = if ile .readline 0 
while line [0:1] != "*": 

nset = line . split [0] 
if nset not in BCnodes : 

BCnodes . append (nset) 
line = if ile .readline 0 

line = if ile .readline 0 

while line and line [0 : 6] .upper () != "*B0UND": 
line = if ile .readline 0 

if ile . closeO 
return BCnodes 

# ******************** wr iteDef ormedInputFile SUB-FUNCTION ******************** 
def writeF ile (copyFile , currInpFile , def _coords , g_def _coords , stepNum, frame , \ 

BCnodes , res) : 

if ile = open(copyFile , "r") # Read the original file for copying 

of ile = open(currInpFile , "w") # Overwrite the original as deformed 

print "***Writing deformed .inp file...\n" 
line = if ile. readline 0 
while line [0 : 5] .upper 0 != "*N0DE": 
of ile .write (line) 
line = if ile .readline 0 

print "len of def_coords:" 

print len(def _coords) 

print "max. val in def_coords:" 

print max(def _coords) 

print "len of g_def _coords : " 

print len(g_def _coords) 

# Write new, deformed nodal information 
of ile .write ("*Node\n") 

for nid in def_coords: 

of ile.write(str(nid)+" , "+str (def _coords [nid] [0] )+" , "+\ 

str (def _coords [nid] [1] )+" , "+\ 
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str (def _coords [nid] [2])+"\n") 


# If local-global analysis, insert the deformed node coords for 

# global region 

line = if ile . readline 0 
if len(g_def _coords) > 0: 

g_currNodeList = [] # A list of global nids to check against 

offset = max(def _coords) # The max. nid of the local region 
nid = line . split [0] 
while nid != offset: 

line = if ile .readline 0 
nid = int (line . split [0] ) 

# After reaching offset, start writing global info 
for i in range (1 ,len(g_def _coords)+l) : 

line = if ile .readline 0 
g_nid = line . split (",") [0] 
g_currNodeList . append (g_nid) 

of ile . write(str (g_nid)+" , "+str (g_def _coords [i] [0] )+" , "+\ 
str (g_def _coords [i] [1] )+" , "+\ 
str (g_def _coords [i] [2])+"\n") 

# Check that the number of global nodes in previous and 

# current models is the same 

if len(g_def _coords) != len(g_currNodeList) : 

print "The number of global nodes do not correspond \ 
between old and new mesh models" 
print "Check script!" 
exit (0) 

# Now copy the rest of the file data 
line = if ile . readline 0 

while line [0 : 8] .upper 0 != "*ELEMENT": 
line = if ile .readline 0 

of ile .write (line) 
line = if ile . readline 0 
while line [0 : 5] .upper 0 != "*STEP": 
of ile. write (line) 
line = if ile .readline 0 

# Write map step before first step 

if res == "Y" : # Step number includes the Equilibrate step \ 

# from first analysis 

newline = "*MAP SOLUTI ON, STEP=7od, INC=7od, UNBALANCED STRESS=RAMP\n**\n" 7o\ 
( St epNum+1 , frame) 
else : 

newline = "*MAP SOLUTI ON, STEP=7od, INC=7od, UNBALANCED STRESS=RAMP\n**\n" 7o\ 
( St epNum, frame) 
of ile . write (newline) 

newline = "** Name: Constrained Type: Displacement/Rotation\n**\n" 
of ile . write (newline) 

newline = "**\n*Step, name=Equilibrate\n*Static\nO . 25 , 1., le-05, 0.25\n" 
of ile . write (newline) 

newline = "**\n*Restart , write, frequency=l, overlay\n" 
of ile . write (newline) 
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newline = "*File Format, ASCII\n" 
of ile . write (newline) 

newline = "*Node File, Frequency=9999\n" 
of ile . write (newline) 
newline = "U,\n" 
of ile . write (newline) 

newline = "**\n** BOUNDARY C0NDITI0NS\n**\n" 

of ile . write (newline) 

for i in range (len (BCnodes) ) : 

newline = "**Name: Fixed_7od Type: Displacement/Rotation\n" % i 
of ile . write (newline) 

newline = "*Boundary\n7oS , 1, l\n7oS, 2, 2\n7oS, 3, 3\n" % \ 

(BCnodes [i] , BCnodes [i] ,BCnodes[i]) 
of ile . write (newline) 

newline = "**\n*End Step\n** \n" 

of ile . write (newline) 

# Write the rest of the file, adding "op=NEW" to boundary 

# condition line, which deactivates BC’s from the mapping step, 
while line: 

if line [0 : 6] .upper 0 == "*B0UND": 
line = "*Boundary, op=NEW\n" 
of ile . write(line) 
line = if ile .readline 0 
continue 

of ile .write (line) 
line = if ile .readline 0 

if ile . closeO 
of ile . close 0 

# * * * * * * * * * * * * * * * * * * * * * * writeDef ormedInputFile MAIN ************************** 

# Move undeformed input file to new directory 

dir = createSubDirectory(currInpFile, "Undeformed") 
fname = currInpFile . split ("/") [-1] . split (". ") 
copyFile = dir+fname [-2] +"_Undef ormed. "+fname [-1] 
if not os .path. exists(copyFile) : 

os . rename (currInpFile , copyFile) 

# Prescan the original input to get nodes (or nodesets) with 

# prescribed BC’s before writing the deformed input file. 

BCnodes = getBCnodesets (copyFile) 

# Write the deformed input file 

writeFile(copyFile , currInpFile ,def_coords ,g_def_coords , stepNum, frame , BCnodes , res) 

if name == " main ": 

print doc 

path = os.getcwd 

res = sys . argv[-10] # "Y" if mapping from a previous restart file 
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localglobal = sys.argv[-9] 
previnp = sys.argv[-8] 
fdbPrev = sys.argv[-7] 
odbFile = sys.argv[-6] 
currInpFile = sys.argv[-5] 
currInpFile_full = sys.argv[-4] 
fdbCurr = sys.argv[-3] 
stepName = sys.argv[-2] 
frame = int (sys . argv [-1] ) 

start_time = time. time () 

# Check for file with initial displacements; get them. 
initial_disps = getInitialDisplacements (f dbPrev) 

# Create MeshTools object of the current and previous . inp files 
prevMesh = createMeshToolsObject (previnp , "INP") 

currMesh = createMeshToolsObject (currInpFile , "INP") 

# Get list of nodes from previous mesh 
oldNodeList = prevMesh. GetNodeList () 

# Function returns stepNum corresponding to stepName, 

# list of all nodes in previous odb, 

# dictionary of odb disps for all nodes **if not local-global** 

# dictionaries of odb disps for local and global regions **if local-global** 

# (empty diets if inapplicable) 

(stepNum,prevFullNodeList ,relative_disps ,l_relative_disps , \ 
g_relative_disps , g_odbIndices) = getRelativeDisplacements (odbFile , \ 

StepName , frame , oldNodeList , localglobal) 


# Previous .inp file (local only, if applicable) 

# local region only if local-global analysis 

# *_full.inp; Empty if not local-global analysis 


if localglobal == "N" : 

# All arguments --> entire mesh 

abs_disps = computeAbsoluteDisplacements (initial_disps , \ 
relative_disps ,prevFullNodeList , prevMesh) 
elif localglobal == "Y" : 

# All arguments --> local region only 

abs_disps = computeAbsoluteDisplacements (initial_disps , l_relative_disps , \ 

oldNodeList , prevMesh) 

# Map displacements (only on local region, if applicable) 

mapped_disps = mapDisplacementsToUndef ormedMesh (prevMesh, abs_disps , currMesh) 
writeMappedDispsToFile (mapped_disps , currInpFile , stepName , frame) 

# Compute deformed coordinates 

def_coords = computeDeformedCoords (currInpFile ,mapped_disps) 
if localglobal == "Y" : 

# Simply add global node disps to node coords of previous odb 
g_def_coords = computeGlobalCoords (odbFile , stepName , frame , \ 

g_relative_disps ,g_odbIndices) 

elif localglobal == "N" : 
g_def_coords = [] 

# Write the deformed input file 
if localglobal == "Y" : 

# Write deformed *_full.inp file; move undeformed input 
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# files to *_Undef ormed folder 

writeDef ormedInpFile(def _coords ,g_def _coords , currInpFile_full , \ 

stepNum, frame , res) 

new_path = str() 
tmp = currInpFile . split ("/") 
fileName = tmp[-l] 
tmp[-l] = "Undeformed" 
for i in range (len (tmp) ) : 
new_path+=tmp [i] +"/" 
new_path+=f ileName 
os . rename (currInpFile ,new_path) 
elif localglobal == "N" : 

# Write deformed * . inp file 

writeDef ormedInpFile(def _coords ,g_def_coords , currInpFile , \ 

StepNum, frame , res) 

printTimeStamp(start_time) 

#exit (0) 


C. Details from Integrally-stiffened Panel Test Program 

Researchers at Alcoa Technical Center fabricated several ISPs as part of a test program in 1998. 
The purpose of the test program was to compare fatigue crack growth and residual strength among 
ISPs machined from either of two, lower wing skin, aluminum alloys — AA 2024-T351 or C433-T39. 
Test information has been utilized by others for analysis purposes [15, 57, 59]. Alcoa Technical 
Center has provided the current authors with test details, which are provided here for completeness. 

Fatigue crack growth and residual strength tests were conducted at Alcoa Technical Center for 
four ISPs, two machined from AA 2024-T351 and two machined from C433-T39. Dimensions of all 
panels are shown in Fig. 16. In each panel, an initial two-bay saw cut of length 2agaw ~ 2.54 cm was 
made at mid-height to completely sever the middle stiffener. The initial cut was then pre-cracked 
to length 2ao ~ 5.08 cm by applying uniaxial cyclic load in the y-direction, where Pmin = 31 kN 
and Pmax = 311 kN. Subsequently, a modified transport wing standard (TWIST) loading spectrum 
[60], shown in Table 6, was applied in the y-direction to propagate the fatigue crack. The applied 
spectrum had a mean flight stress of Smf = 08.9 MPa, truncated to level V, and included a ground- 
air-ground (GAG) cycle with a reduced ground level stress of Sground = —34.5 MPa. Taxi loads 
were neglected. Depending on the panel, the fatigue crack was propagated until both crack fronts 
were 2.54 cm short of reaching the intact stiffeners (2a^ ~ 24.1 cm) or until both crack fronts just 
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entered the intact stiffeners {2ai ^ 30.0 cm). Subsequently, each panel was loaded monotonically in 
uniaxial tension until failure occurred by unstable crack growth. A test matrix of the four panels 
and corresponding residual strengths is shown in Table 7. The ISP simulated in this work follows 
panels 5 and 5 A. 


Table 6 Standard TWIST spectrum scaled to mean flight stress Smf — 68.9 MPa with variable 
amplitude loads Sa (left). Modified spectrum (right) applied uniaxially to the Alcoa ISP prior 
to residual strength testing. Courtesy Alcoa Technical Center. 

Standard TWIST profile scaled to Modified TWIST profile applied to 
S^f = 68.9 Mpa ISPs 


Block 

Level 

Load Level 

Number of 

(Smf ± Sa ), 

cycles per 

MPa 

4000 flights 

I 

68.9 ± 110.3 

1 

II 

68.9 ± 103.4 

2 

III 

68.9 ± 89.6 

5 

IV 

68.9 ± 79.3 

18 

V 

68.9 ± 68.6 

52 

VI 

68.9 ± 57.9 

152 

VII 

68.9 ±47.2 

800 

VIII 

68.9 ± 36.5 

4170 

IX 

68.9 ± 25.9 

34800 

X 

68.9 ± 15.3 

358665 

GAG 

-34.5 (last 
peak per flight) 

4000 


Block 

Level 

Load Level 
(Smf ± Sa ), 
MPa 

Number of 
cycles per 
4000 flights 

I 

- 

truncate loads 
to level V 

II 

- 

III 

- 

IV 

- 

V 

68.9 ± 68.6 

78 

VI 

68.9 ± 57.9 

152 

VII 

68.9 ±47.2 

800 

VIII 

68.9 ± 36.5 

4170 

IX 

68.9 ± 25.9 

34800 

X 

68.9 ± 15.3 

358665 

GAG 

-34.5 (last 
peak per flight) 

4000 


Table 7 Description of all ISPs and corresponding residual strengths from Alcoa test program. 
Courtesy Alcoa Technical Center. 


Panel ID 

Material 

Crack length, 2a,, 
just prior to residual 
strength test (cm) 

Residual 
strength (kN) 

3 

2024-T351 

24.1 cm 

1225 

4 

2024-T352 

30.0 cm 

1015 

5A 

C433-T39 

24.1 cm 

1660 

5 

C433-T39 

24.1 cm 

1385 

6 

C433-T40 

30.0 cm 

1295 
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